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1  INTRODUCTION 


The  subject  of  this  report  is  traveling  waves  with  complex  propagation  constants  (complex  waves) 
on  one-,  two-,  and  three-dimensional  (ID,  2D,  3D)  periodic  arrays  of  lossless  and  lossy  magnetodi¬ 
electric  spheres.1  Such  arrays  are  an  important  class  of  artificial  composite  materials  materials 
formed  by  embedding  inclusions  of  one  material  in  a  matrix  of  another  kind  of  material.  Just  like 
natural  materials,  artificial  composite  materials  can  have  properties  that  are  very  different  from 
those  of  their  constituents  when  the  scale  of  use  or  observation  of  the  composite  material  (the 
macro-level)  is  much  larger  than  the  scale  of  the  detailed  structure  of  the  composite  material  (the 
micro-level).2 3  For  example,  it  has  been  known  since  the  seventeenth  century  that  ruby-colored  glass 
can  be  made  by  adding  nanometer-sized  gold  particles  to  silicate  glass  [1].  Since  the  1940's  much 
attention  has  been  given  to  artificial  dielectrics  made  from  lattices  of  different  kinds  of  conducting 
elements.  Similar  interest  in  artificial  magnetic  materials  dates  from  the  1950's.1  The  subject 
of  photonic  crystals  [6]-[8]  is  a  recent  and  very  important  extension  of  earlier  work  on  artificial 
dielectrics  to  the  optics  frequency  range.  More  recently  there  has  been  great  interest  in  artificial 
composite  materials  characterized  by  a  negative  index  of  refraction  (group  and  phase  velocities  in 
opposite  directions),  especially  those  having  both  their  bulk  permittivity  and  permeability  nega¬ 
tive  [double  negative  (DNG)  materials]  [9]-[18].4  In  addition  to  the  term  “metamaterials”  being 
commonly  applied  to  such  artificial  negative  index  or  doubly  negative  materials,  the  term  has  also 
come  to  be  applied  to  any  materials  not  found  in  nature,  such  as  artificial  dielectric  and  magnetic 
materials,  and  photonic  crystals,  whose  microstructures  have  been  engineered  to  result  in  desired 
macro-properties.  The  most  interesting  macro-  or  bulk  properties  of  electromagnetic  metamaterials 
are  those  related  to  how  electromagnetic  waves  interact  (scatter,  propagate,  refract,  etc.)  with  the 
materials,  hence  the  central  importance  in  this  report  of  the  waves  that  can  be  supported  by  the 
periodic  arrays  we  study.  Indeed,  as  will  be  shown  below,  when  the  free-space  and  traveling-wave 
electrical  separation  distances  between  the  elements  in  a  3D  cubic-lattice  periodic  array  are  both 
small,  the  array  can  be  regarded  as  a  homogeneous  isotropic  medium  whose  bulk  permittivity  and 
permeability  can  be  obtained  from  the  solution  to  the  dispersion  equation  for  the  array.  (See, 
however,  Footnote  23  in  Section  5.) 

Our  investigation  of  complex  traveling  waves  on  periodic  arrays  of  magnetodielectric  spheres 
is  motivated  in  part  by  the  theoretical  demonstration  by  Holloway  et  al.  [23],  based  on  work  by 
Lewin  [24],  that  a  DNG  material  can  be  formed  by  embedding  an  array  of  spherical  magnetodieletric 

]The  term  “leaky  wave”  is  sometimes  used  synonymously  with  the  term  “complex  wave.”  Nonetheless,  in  this 
report  we  reserve  the  term  “leaky  wave”  for  outgoing  improper  complex  waves;  see  Appendix  C. 

2  Almost  needless  to  say,  the  distinction  between  macro  properties  and  microstructure  is  so  basic  and  far  reaching 
as  to  be  taken  for  granted  in  our  world  view  of  the  molecular  and  atomic  structure  of  matter  or  of  the  cellular 
structure  of  biological  organisms. 

3 A  good  summary  of  early  work  in  artificial  dielectrics  and  magnetics  is  given  in  [2].  A  review  of  early  work  in 
artificial  plasmas  as  well  as  of  early  work  in  artificial  dielectrics  and  magnetics  can  be  found  in  [3].  A  static-field 
treatment  of  artificial  dielectrics  is  given  in  [4]  and  [5,  ch.  12]. 

4 The  chapter  in  [17]  by  Simovski  and  Tretyakov  referenced  above  (see  [3])  includes  a  historical  overview  of  early 
work  on  backward  waves  and  negative  refraction,  and  an  excellent  discussion  of  the  controversy,  dating  from  the 
early  2000’s,  regarding  the  possibility  of  negative  refraction.  Simovski  and  Tretyakov  conclude  that  all  the  arguments 
that  have  been  advanced  to  disprove  the  possibility  of  negative  refraction  have  been  convincingly  refuted.  (In  this 
connection,  it  is  unfortunate  that  the  recent  posthumous  publication  of  a  book  by  Munk  [19]  has  revived  arguments 
against  the  possibility  of  negative  refraction  and  DNG  materials  that  were  refuted  years  ago.  Indeed  our  own  work  [20, 
sec.  12. 5], [21,  sec.  6.4], [22]  provides  clear  examples  of  theoretically  possible  DNG  metamaterials.  While  skepticism 
regarding  the  practical  importance  of  DNG  structures  in  the  foreseeable  future  is  fully  justified  considering  the  paucity 
of  significant  practical  applications  to  date,  the  validity  of  the  theoretical  foundations  of  negative  refraction  and  DNG 
metamaterials  is  well  established.) 
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particles  in  a  background  matrix5  and  by  the  current  widespread  interest  in  DNG  metamaterials. 
An  important  application  of  our  work  is  to  the  design  of  composite  DNG  metamaterials.  Our  work 
can  be  used,  for  example,  to  investigate  how  the  frequency  regions  in  which  an  array  behaves  as 
an  isotropic  homogeneous  DNG  medium  depend  on  the  permittivity,  permeability,  and  packing 
of  the  array  elements,  and  to  estimate  how  losses  in  the  elements  affect  the  attenuation  of  waves 
propagating  in  composite  DNG  structures.  The  applicability  of  our  work  is  not  restricted  to 
composite  structures,  however.  It  is  also  relevant,  for  example,  to  understanding  the  propagation 
of  plasmonic  waves  on  periodic  chains  of  metallic  optical  nanospheres,  to  linear  arrays  of  electric 
or  magnetic  dipoles  such  as  the  Yagi-Uda  antenna,  and  to  modeling  frequency  selective  surfaces 
composed  of  planar  periodic  arrays  of  small  electrically  or  magnetically  polarizable  elements. 

The  work  to  be  described  here  is  an  extension  of  our  previous  investigations  of  traveling  waves 
with  real  propagation  constants  supported  by  ID,  2D,  and  3D  periodic  rectangular-lattice  arrays 
of  lossless  acoustic  monopoles,  electric  or  magnetic  dipoles,  and  magnetodielectric  spheres  [25]- 
[30], [20]- [22]  using  a  spherical- wave  source  scattering-matrix  formulation.  Unlike  our  previous  work, 
in  this  report  the  array  elements  can  be  either  lossless  or  lossy,  and  the  propagation  constants  can 
be  either  real  or  complex.  The  propagation  constants  are  always  complex  when  the  array  elements 
are  lossy,  and  can  also  be  complex  even  when  the  array  elements  are  lossless,  as  in  the  fast-wave 
region  of  ID  and  2D  arrays  of  lossless  magnetodielectric  spheres,  occasionally  in  the  slow- wave 
region  of  these  arrays  close  to  the  light  line,6  and  in  the  stop-band  regions  (where  no  traveling 
wave  with  real  (3d  exists)  of  ID,  2D,  and  3D  arrays.  Additionally,  complex  eigenmodes  with  high 
attenuation  can  be  supported  by  ID,  2D,  and  3D  periodic  arrays  of  lossless  scatterers. 

The  focus  of  our  attention  is  the  so-called  k-(3  equation  (diagram)  —  in  our  work  more  properly 
referred  to  as  the  kd-(3d  equation  (diagram)  that  relates  the  traveling  wave  electrical  separation 
distance  (3d  of  the  array  elements  in  the  direction  parallel  to  a  specified  array  axis  along  one  of 
the  principal  rectangular-lattice  directions,  to  the  corresponding  free-space  electrical  separation 
distance  kd,  where  k  =  uj/c  is  the  free-space  wavenumber  with  uj  >  0  the  angular  frequency  and  c 
the  free-space  speed  of  light.  An  implicit  harmonic  time  dependence  of  exp(-iu;f)  is  assumed  here 
and  throughout  the  report. 

In  this  report  we  treat  traveling  waves  only  in  the  direction  parallel  to  a  specified  array  axis. 
Our  objective  is  to  formulate  an  efficient  rigorous  method  to  obtain  the  traveling- wave,  array- 
element  electrical  separation  distance  /3d  (“propagation  constant”)  of  all  the  discrete  waves  -  real 
or  complex,  proper  or  improper  -  propagating  parallel  to  the  array  axis  that  can  be  supported  by 
the  array  for  a  particular  free-space,  array-element  electrical  separation  distance  kd  (“frequency”). 
This  is  accomplished  by  deriving  and  solving  the  exact  (to  within  the  dipole  approximation)  source- 
free  kd  [3d  equation  for  the  array.  Thus,  all  the  “propagation  constants”  of  the  discrete  waves  are 
found  from  an  analysis  of  the  source-free  array,  without  having  to  resort  to  determining  the  full 

5Lewin’s  procedure  for  obtaining  the  effective  permittivity  and  permeability  of  a  mixture  consisting  of  a  homoge¬ 
neous  material  in  which  is  embedded  a  cubic  lattice  of  magnetodielectric  particles  is  quite  different  from  ours  although 
there  are  some  similarities.  He  considers  the  problem  of  a  plane  wave  incident  on  a  half-space  of  the  mixture,  and 
obtains  the  propagation  constant  of  the  transmitted  wave  by  first  calculating  the  electric  and  magnetic  field  incident 
on  a  reference  particle  as  the  sum  of  the  transmitted  wave  plus  the  sum  of  the  Mie  electric  and  magnetic  dipole  fields 
scattered  from  all  the  particles  of  the  lattice  other  than  the  reference  particle  itself.  He  then  evaluates  the  summations 
by  approximating  them  by  integrals,  leading  to  a  simultaneous  pair  of  integral  equations  which  can  be  solved  for 
the  propagation  constant  of  the  transmitted  wave.  This  in  turn  allows  him  to  calculate  the  wave  reflected  from  the 
interface,  and  knowing  both  the  transmitted  and  reflected  waves  he  obtains  expressions  for  the  bulk  parameters  of 
the  mixture.  Although  ingenious,  his  procedure  is  far  more  complicated,  unnecessarily,  than  ours  and  is  approximate 
rather  than  exact  as  ours  is  (to  within  the  dipole  approximation).  Furthermore,  unlike  Lewin’s  work,  in  our  work  the 
conditions  for  the  validity  of  the  bulk  parameter  expressions  we  derive  are  extremely  simple. 

6In  our  plots  of  kd-(3d  (dispersion)  diagrams  the  light  line  is  the  line  3 ft(/3d)  =  kd ,  to  the  left  of  which  is  the 
fast- wave  region,  0  <  $t((3d)  <  kd ,  and  to  the  right  of  which  is  the  slow- wave  region,  kd  <  %t(/3d)  <  tt. 
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Green’s  function  solution  for  the  array.' 

The  restriction  to  waves  traveling  in  a  direction  parallel  to  a  specified  array  axis  greatly  sim¬ 
plifies  the  analysis  because  of  the  rectangular-lattice  periodic  arrays  investigated  here;  the  analysis 
of  waves  traveling  in  an  arbitrary  direction  in  2D  or  3D  arrays  would  be  far  more  complicated. 
However,  even  though  the  traveling  waves  treated  are  assumed  to  be  in  a  direction  parallel  to  an 
array  axis,  it  will  be  seen  below  (see  Section  5)  that  when  both  \/3d\  and  kd  are  less  than  about  one, 
a  3D  cubic-lattice  array  of  identical  magnetodielectric  spheres  can  be  regarded  as  a  homogeneous 
isotropic  medium  whose  bulk  permittivity  and  permeability  are  obtainable  from  the  solution  of  the 
kd-(3d  equation  for  a  transverse  wave  traveling  parallel  to  the  array  axis.  Thus,  the  solution  of  the 
dispersion  equation  of  a  wave  propagating  parallel  to  the  direction  of  an  axis  of  a  3D  cubic-lattice 
array  leads  directly  to  knowing  in  what  frequency  region  the  array  can  be  regarded  as  a  homoge¬ 
neous  isostropic  medium,  and  to  knowing  the  values  of  the  bulk  parameters  in  that  frequency  region 
determining  wave  propagation  in  an  arbitrary  direction  in  that  medium.  (See,  however,  Footnote 
23  in  Section  5.) 

The  class  of  problems  analyzed  in  this  report  can  be  described  as  follows.  We  consider  periodic 
arrays  -  1 D,  2D,  or  3D  of  identical,  radially  symmetric  elements,  lossy  or  lossless,  with  the  elements 
characterized  by  two  scattering  coefficients  that  relate  the  electric  and  magnetic  elementary  dipole 
fields  scattered  from  an  element  to  the  incident  electric  and  magnetic  fields  at  the  element  center. 
As  in  our  previous  related  work,  it  is  assumed  that  the  array  elements  are  sufficiently  small,  or  the 
frequency  ranges  investigated  sufficiently  low,  so  that  only  the  fields  of  the  lowest  order  spherical 
multipoles  (electromagnetic  dipoles)  are  significant  in  analyzing  scattering  from  the  array  elements. 
All  the  multiple  interactions  among  these  dipole  elements  are  taken  into  account  so  that  our  analyses 
are  exact  to  within  the  dipole  field  approximation.  The  spacing  of  elements  in  the  direction  parallel 
to  the  array  axis  is  denoted  by  d  and  in  the  direction  or  directions  normal  to  the  array  axis  by  h.  For 
any  of  these  arrays  we  obtain  and  solve  the  kd  fid  equation,  where  (3d  can  be  either  real  or  complex, 
and  plot  the  resulting  kd-(3d  diagram.  A  3D  array,  when  both  \(3d\  and  kd  are  less  than  about 
one,  and  when  the  longitudinal  element  separation  d  equals  the  transverse  element  separation  h , 
can  generally  be  regarded  as  a  homogeneous  isotropic  medium  with  effective  bulk  permittivity  and 
permeability.  When  |/?d|  and  kd  are  sufficiently  small  for  a  3D  array  to  be  regarded  as  homogeneous 
and  isotropic,  we  obtain  values  for  these  bulk  parameters  from  the  solution  of  the  kd  (id  equation 
for  the  array.  (See,  however,  Section  5,  Footnote  23.) 

Some  basic  properties  of  the  kd-(3d  diagram  may  be  noted  here.  In  [35,  ch.  7]  it  is  shown  that 
the  dependence  of  the  kd~(3d  diagram  on  (id  is  periodic  in  (3d  with  a  period  of  2n.  In  Appendix 
A  of  [20]  and  the  paper  [36],  it  is  proven  that  if  a  periodic  array  of  reciprocal  elements  supports  a 
traveling  wave  with  propagation  constant  (i,  it  also  supports  a  corresponding  traveling  wave  with 
propagation  constant  -(3.  Therefore,  for  periodic  arrays  of  reciprocal  elements,  as  are  all  the  arrays 
considered  in  this  report,  kd  is  an  even  function  of  (3d.  It  follows  that  for  ?R{(3d)  in  the  interval 
7 r  <  ft{(3d)  <  2tt 

kd((3d)  =  kd{ 2tt  -  (id),  7 r  <  5?(/3d)  <  2tt  (1.1) 

7 The  restriction  of  our  analysis  to  source- free  arrays  has  the  limitation  that  it  cannot  be  used  to  solve  problems 
such  as  those  treated,  for  example,  in  [31]-[34],  involving  periodic  arrays  excited  by  sources. 
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where  we  have  written  kd  as  a  function  of  (3d.  Hence  we  need  consider  / 3d  only  in  the  interval* 

0  <  $R(/3d)  <  tt  .  (1.2) 

If  we  restrict  5R(/?d)  to  this  positive  interval,  however,  G(/?d)  needs  to  range  over  negative  as  well 
as  positive  values  in  order  to  find  all  the  discrete  complex  waves  that  can  propagate  on  arrays  with 
lossy  elements;  see  Appendix  C. 

Although  there  is  no  upper  limit  on  the  transverse  inter-element  spacing  h  for  either  2D  or 
3D  periodic  arrays,  the  expressions  we  give  for  the  rapidly  convergent  summations  in  the  kd-(3d 
equations  are  valid  only  for  kh  <  2i r,  that  is,  for  h  less  than  a  wavelength.  This  restriction  on 
the  size  of  h  is  not  an  essential  limitation  of  either  the  transverse  element  separation  or  of  the 
analyses  we  perform.  It  is,  rather,  a  matter  of  our  not  wanting  to  unnecessarily  complicate  the 
form  of  the  rapidly  convergent  expressions  we  give  by  making  them  independent  of  the  range  of 
kh  since  in  most  practical  applications  the  transverse  element  spacing  can  be  expected  to  be  less 
than  a  wavelength.  Examples  of  how  the  range  of  kh  can  be  extended  can  be  found  in  [20,  secs.  2 
and  3].  Also  of  interest  are  the  limiting  values  of  the  kd-(3d  equations  as  kh  — >  2tt  since  some  of 
the  individual  terms  of  the  rapidly  convergent  expressions  in  the  kd  (3d  equations  are  singular  at 
kh  =  2n.  Closer  analysis  shows,  however,  that  the  singularities  of  the  various  terms  in  the  kd  (3d 
equation  cancel  one  another  and  hence  the  kd-(3d  equations  remain  non-singular  at  kh  =  2n . 

We  follow  the  same  procedure  for  all  the  arrays  considered  in  this  report.  We  begin  as  in 
[27], [29], [20], [21]  by  assuming  the  existence  of  a  wave  with  real  propagation  constant  (3  traveling 
from  —  oo  to  -hoc  in  the  direction  parallel  to  the  array  axis.  We  are  not  concerned  in  any  way 
(other  than  in  Appendix  B)  with  the  excitation  of  the  wave.  No  explicit  source  for  this  traveling 
wave  is  assumed  (even  though,  of  course,  the  wave  cannot  exist  without  some  source;  see  Appendix 
B).  We  simply  suppose  that  such  a  traveling  wave  supported  by  the  array  elements  exists  and, 
given  that,  proceed  as  follows  to  find  the  kd  (3d  equation  expressing  (3d  as  a  function  of  kd  and 
the  parameters  (permittivity,  permeability,  and  packing)  of  the  array  elements.  The  traveling 
wave  excites  the  array  elements  which  then  radiate  scattered  dipole  fields  which  in  turn  also  serve 
as  excitations  for  the  array  elements.  (More  precisely,  the  traveling  wave  is  the  resultant  of  the 
multiplicity  of  spherical  vector  dipole  waves.)  Orthogonality  of  the  vector  spherical  wave  functions 
and  the  radially  symmetric  scatterers  assumed,  imply  that  an  incident  electric  (magnetic)  dipole 
field  can  excite  only  a  scattered  electric  (magnetic)  dipole  field  [27,  p.  6].8 9  Additionally,  noting 
that  the  vector  spherical  dipole  fields  are  the  only  multipole  vector  spherical  fields  that  are  non¬ 
zero  at  an  element  center,  an  incident  electric  (magnetic)  field  at  an  element  center  results  in  a 
scattered  electric  (magnetic)  dipole  field  radiated  by  an  elementary  electric  (magnetic)  dipole  in  the 
direction  of  the  incident  electric  (magnetic)  field.  The  coefficient  of  the  scattered  electric  (magnetic) 
dipole  field  is  given  by  the  product  of  the  complex  electric  (magnetic)  dipole  scattering  coefficient 
and  the  coefficient  of  the  incident  electric  (magnetic)  dipole  field.  When  the  array  elements  are 
homogeneous  spheres  the  dipole  scattering  coefficients  are  the  ordinary  Mie  electric  and  magnetic 
dipole  coefficients,  given,  for  example,  by  [38,  sec.  9.25  (10), (11)].  More  complicated  elements  such 

8 For  leaky  wave  antennas  (LWA’s)  formed  by  periodically  placed  holes  (perforations)  in  the  side  of  an  otherwise 
uniform  closed  waveguide,  some  authors  [37]  define  the  real  propagation  constant  po  with  n  <  pod  <  in  the 
imperforated  waveguide  as  the  primary  (n  =  0)  Floquet-mode  propagation  constant  (even  though  it  does  not  describe 
the  variation  of  the  periodic  fields  outside  the  perforated  waveguide).  In  that  case,  our  %t(pd)  corresponds  to  their 
pod  —  2n.  That  is,  their  n  =  —  1  space  harmonic  for  such  LWA’s  equals  our  n  =  0  space  harmonic  (our  primary 
Floquet  mode). 

®The  vector  spherical  wave  functions  underlying  our  analysis  are  given  by  [29,  eqs.  (l)-(3)],[27,  eqs.  (16a,b), 
eqs.  (64), (65)].  In  this  representation,  an  electric  (magnetic)  dipole  wave  incoming  on  a  sphere  has  only  an  electric 
(magnetic)  field  at  the  sphere  center. 
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as  coated  spheres  require  appropriate  scattering  coefficients.10  We  sum  the  electromagnetic  fields 
incident  on  a  reference  element  from  all  the  other  elements  of  the  array,  the  summations  taking 
into  account  all  the  multiple  scattering  interactions  of  the  array  elements.  The  field  incident  on 
a  reference  element  from  all  the  other  elements  of  the  array  is  called  the  “interaction  field”  in 
the  literature  [5,  ch.  12], [39,  sec.  4. 5.1].*  11  In  the  case  of  2D  and  3D  arrays  we  then  replace  the 
summations  over  the  indices  of  the  axes  transverse  to  the  array  axis  by  expressions  obtained  in  [20] 
to  somewhat  increase  the  rate  of  convergence  of  the  original  summations.  At  this  point  we  have  in 
effect  an  equation,  containing  infinite  summations  over  the  index  of  the  array  axis,  that  determines 
fid  as  a  function  of  kd  and  the  array  element  parameters. 

When  (id  is  real,  these  infinite  summations  converge,  but  very  slowly.  However,  rapidly  conver¬ 
gent  expressions  can  be  found  for  them  yielding  the  kd  [id  equation  in  a  form  that  can  be  efficiently 
evaluated  and  solved.  If  fid  is  complex,  however,  then  this  final  step  in  the  procedure  for  obtaining 
the  kd-fid  equation  for  traveling  waves  with  real  propagation  constants  supported  by  periodic  arrays 
of  lossless  scatterers  can  no  longer  be  taken  because  the  summations  over  the  index  of  the  array  axis 
diverge.  (This  apparent  paradox  of  the  divergence  of  these  summations  from  which  the  complex 
propagation  constants  are  obtained  indirectly  is  discussed  in  Appendix  B.)  Accordingly,  following 
Citrin  [41],  Koenderink  [42],  and  Alii  and  Engheta  [43],  we  rewrite  the  kd-fid  equation  (containing 
the  divergent  summations)  in  a  form  that  can  be  analytically  continued  into  the  complex  fid  plane 
and  evaluated  for  any  fid.12  This  equation  reduces  to  the  kd-fid  equation  obtained  in  our  investi¬ 
gations  of  traveling  waves  with  real  propagation  constants  on  lossless  arrays  when  (id  is  real,  and 

,0Because  of  the  symmetry  of  the  array  lattice,  for  traveling  waves  with  the  electric  and  magnetic  fields  perpendic¬ 
ular  to  the  traveling- wave  propagation  direction  (the  array  axis),  all  components  of  the  electric  and  magnetic  fields 
incident  on  the  center  of  an  element  cancel  except  for  an  electric  and  magnetic  field  perpendicular  to  each  other  and 
to  the  array  axis,  and  hence  the  array  elements  effectively  behave  as  pairs  of  crossed  electric  and  magnetic  dipoles 
perpendicular  to  each  other  and  to  the  direction  of  propagation  of  the  traveling  wave.  This  should  not  be  understood, 
however,  as  meaning  that  the  array  elements  are  in  general  pairs  of  crossed  electric  and  magnetic  dipoles  that  are 
fixed  independently  of  the  dipole  fields  incident  on  the  elements.  In  this  report  we  will  refer  to  “effective”  dipoles  to 
emphasize  this.  Since  our  analysis  is  for  waves  traveling  in  the  direction  of  an  axis  of  a  rectangular  lattice  array,  it  can 
apply  not  only  to  radially  symmetric  elements  with  no  fixed  electric  and  magnetic  dipoles,  but  also  to  elements  that 
can  only  be  polarized  electrically  and  magnetically  in  fixed  directions  perpendicular  to  each  other.  Indeed,  we  can 
obtain  all  our  results  by  replacing  the  radially  symmetric  elements  with  no  fixed  polarization  direction  by  elements 
with  such  fixed  polarizability.  But  physically  the  two  kinds  of  elements  are  not  at  all  identical.  If  the  electric  and 
magnetic  polarizability  directions  are  fixed  independently  of  the  incident  dipole  fields,  then  for  our  analysis  to  be 
applicable  the  induced  electric  and  magnetic  dipoles  must  be  uncoupled  so  that  an  incident  electric  (magnetic)  field 
at  the  element  center  excites  only  the  electric  (magnetic)  dipole  field. 

1 1  For  1-D  and  2-D  arrays,  the  summations  giving  the  interaction  field  for  real  (id  converge,  the  2-D  array  summations 
more  slowly  than  the  1-D  summations.  For  3-D  arrays,  however,  it  is  necessary  to  add  a  small  imaginary  part  to 
fid  if  the  summations  are  to  converge  rather  than  oscillate.  Additionally,  the  3-D  interaction  field  can  be  viewed 
as  a  summation  of  the  contributions  from  the  planes  of  elements  perpendicular  to  the  array  axis.  The  summations 
over  the  transverse  direction  indices  to  obtain  the  contributions  to  the  interaction  field  of  these  planes  of  elements 
converge,  but  extremely  slowly;  see  the  integration  analogue  in  [40,  sec.  D.2  of  Appendix  D],  In  any  event,  for  all 
arrays,  1-D,  2-D,  and  3-D,  direct  summation  of  the  scattered  dipole  fields  is  completely  impractical  for  obtaining 
kd-fid  diagrams,  and  hence  the  conversion  of  the  original  summations  to  rapidly  convergent  forms  is  central  to  the 
analysis  of  traveling  waves  supported  by  the  arrays. 

12Complex  propagation  constants  0  on  infinite  uniform  or  periodic  waveguides  arise  from  the  solution  for  the  fields 
that  are  excited  by  sources  illuminating  these  infinite  waveguides.  Specifically,  the  fields  along  the  longitudinal  2  axis 
of  an  infinite  waveguide  can  be  expressed  as  a  Fourier  transform  over  real  values  of  0  ranging  from  -00  to  +00.  If 
this  integration  along  the  infinite  real  (i  line  is  evaluated  by  closing  the  contour  in  the  complex  li  plane,  the  residue 
contributions  from  the  poles  of  the  integrand  represent  waves  with  complex  propagation  constants  0.  The  poles  of 
this  integrand  are  identical  to  the  zeros  of  the  homogeneous  (source-free)  transcendental  k-0  equation  obtained  for 
real  (i  and  analytically  continued  into  the  complex  0  plane;  see  Appendix  A.  Thus,  we  are  able  to  obtain  the  complex 
propagation  constants  for  the  complex  waves  excited  by  sources  illuminating  a  waveguide  from  the  transcendental 
equation  for  the  real  propagation  constants  on  the  source-free  waveguide. 


5 


enables  complex  values  of  fid  to  be  found  as  a  function  of  kd ,  kh ,  and  the  array  element  parameters 
when  the  array  elements  are  lossy,  or  when  we  are  in  the  fast-wave  region,  0  <  3f(/3d)  <  kd,  (and 
even  in  the  slow- wave  region  close  to  the  light  line)  for  ID  and  2D  arrays  of  lossless  elements,  or 
when  we  are  in  the  stop-band  regions  of  ID,  2D,  and  3D  arrays  of  lossless  scatterers.  Additionally, 
complex  eigenmodes  with  high  attenuation  (large  Q(/M))  can  be  supported  by  ID,  2D,  and  3D 
periodic  arrays  of  lossless  scatterers.  For  ID  arrays,  the  analytically  continuable  form  of  the  kd-/3d 
equation  is  obtained  by  replacing  Clausen  functions  by  polylogarithm  functions.  For  2D  arrays,  the 
analytically  continuable  form  of  the  kd  fid  equation  is  found  by  obtaining  analytically  cont  inuable 
and  rapidly  convergent  expressions  for  the  Schlomilch  series13  with  real  fid  that  arise  in  the  course 
of  converting  the  summations  of  the  initial  form  of  the  kd-fid  equation  to  rapidly  convergent  forms. 
This  involves  a  careful  choice  of  branch  cuts.  Obtaining  an  analytically  continuable  form  of  the 
kd  -fid  equation  is  simplest  for  the  3D  arrays  for  which  all  that  is  necessary  is  the  use  of  cos  (fid)  and 
sin(/?d)  when  fid  is  complex.  Much  the  same  conclusion  applies  to  the  3D  alternating  two-sphere 
arrays  treated  in  [22],  and  the  analysis  there  has  been  written  for  complex  waves  and  lossy  arrays 
from  the  start. 

Considerable  insight  into  the  analytic  continuation  procedure,  the  key  step  in  all  the  analysis 
of  this  report  ,  can  be  obtained  from  the  analysis  performed  in  Appendix  A  of  waves  supported  by 
a  2D  slab  of  magnetodielectric  material.  In  that  appendix  it  is  shown  that  analytically  continuing 
the  homogeneous  kd-fid  equation  for  real  fid  into  the  complex  plane,  is  equivalent  to  evaluating 
the  Fourier  transform  expression  for  the  solution  of  the  inhomogeneous  slab  problem  by  closing  the 
Fourier  transform  contour  of  integration  along  the  real  axis  in  the  complex  plane  and  calculating  the 
residues  at  the  zeroes  of  the  integrand.  We  emphasize  that  this  analytic  continuation  procedure 
allows  us  to  find  the  propagation  constants  of  all  the  complex  waves  supported  by  the  arrays, 
improper  as  well  as  proper,  without  any  consideration  of  the  excitation  of  the  waves.  Green’s 
functions  are  not  required  to  find  the  propagation  constants  of  the  complex  waves  supported  by 
the  arrays. 

It  is  important  to  comment  on  the  procedure  used  to  solve  the  kd-fid  equations  for  complex  fid. 
In  our  previous  work  on  traveling  waves  with  real  propagation  constants,  all  the  kd-fid  equations 
encountered  were  real  functions  of  the  real  argument  fid  and  so  a  simple  search  procedure  with 
secant  algorithm  refinement  sufficed.  In  the  kd-fid  equations  encountered  here  in  our  extension  of 
our  earlier  work  to  traveling  waves  with  complex  propagation  constants,  the  equations  are  com¬ 
plex  functions  of  the  complex  argument  fid  and  the  procedure  for  solving  the  equations  for  fid  is 
necessarily  more  complicated.  Starting  with  an  idea  of  a  relatively  large  region  in  the  complex  fid 
plane  where  the  zero  of  the  kd-fid  equation  is  thought  to  lie  (this  can  often  be  obtained  by  starting 
with  the  kd  fid  diagram  for  the  real  case),  we  regard  the  function  of  fid ,  f(fid,  kd ),  whose  zero  we 
are  attempting  to  find,  as  a  function  of  the  two  real  variables  S l{fid),  9( fid )  and  vary  the  real  and 
imaginary  parts  of  fid  separately  with  a  relatively  coarse  grid,  to  locate  the  point  in  the  region  that 
gives  a  minimum  value  of  \f(fid,  kd)  |.  The  search  region  and  grid  are  then  reduced  in  size  and  the 
procedure  repeated  to  obtain  a  more  refined  estimate  of  the  zero.  The  value  ot  fid  that  minimizes 
f(fid,  kd)  obtained  after  several  repetitions  of  this  procedure  is  then  used  as  the  starting  guess  for  a 
quasi-Newton  method  with  a  finite-difference  gradient  [44.  Appendix  A]  implemented  in  the  IMSL 
subroutine  UMINF,  or  for  Muller’s  Method  [45], [46]  implemented  in  the  IMSL  subroutine  ZANLY, 
to  minimize  \f(fid,kd)\  and  thereby  obtain  a  highly  accurate  value  for  the  zero  of  f(fid,kd).  It 
should  be  noted  that  it  is  difficult  to  be  certain  that  all  complex  zeroes  of  the  kd-fid  equation  have 
been  found  in  a  given  region  of  the  complex  fid  plane,  so  that  not  all  kd-fid  curves  may  be  present 
in  the  kd-fid  diagrams  shown  in  Section  8.  A  detailed  description  of  a  procedure  for  finding  the 

13 Series  of  the  form  EajZ„(jx)  where  a3  can  be  a  trigonometric  function. 


complex  modes  for  determining  the  electromagnetic  fields  of  Hertzian  dipoles  in  layered  media  is 
given  in  [47],  but  the  dispersion  equation  for  this  problem  is  considerably  simpler  than  the  equa¬ 
tions  encountered  in  our  analyses  of  periodic  arrays,  and  the  applicability  of  the  procedure  in  [47] 
to  our  problems  is  questionable. 

The  outline  of  the  report  is  as  follows.  There  are  two  main  parts  of  the  report.  The  first 
part,  Sections  2  to  7,  is  devoted  to  the  analysis,  and  the  second  part,  Section  8,  to  the  numerical 
results.  In  Section  2  we  derive  the  kd  (3d  equations  for  the  complex  propagation  constants  on  waves 
supported  by  ID  arrays  of  magnetodielectric  spheres,  lossless  or  lossy.  Subsection  2.1  treats  arrays 
with  the  effective  electric  and  magnetic  dipoles  of  the  spheres  oriented  at  right  angles  to  each  other 
and  to  the  array  axis.  Subsection  2.2  considers  arrays  with  the  dipoles  aligned  with  the  array  axis. 
Since  in  this  case  there  is  no  coupling  between  the  electric  and  magnetic  dipoles,  this  is  equivalent 
to  considering  ID  arrays  of  either  electric  or  magnetic  dipoles  aligned  with  the  array  axis. 

Section  3  is  devoted  to  the  kd  (3d  equations  for  the  complex  propagation  constants  of  waves 
supported  by  2D  arrays  of  lossless  or  lossy  magnetodielectric  spheres.  Subsection  3.1  deals  with 
arrays  with  the  effective  electric  and  magnetic  dipoles  of  the  spheres  at  right  angles  to  each  other 
and  to  the  array  axis.  There  are  two  polarizations  to  be  considered.  In  3.1.1  the  electric  dipoles  are 
taken  to  be  in  the  plane  of  the  array  and  the  magnetic  dipoles  are  perpendicular  to  the  array  plane, 
while  in  3.1.2  the  electric  dipoles  are  perpendicular  to  the  array  plane  and  the  magnetic  dipoles 
are  in  the  plane  of  the  array.  The  kd-(id  equation  for  this  latter  polarization  can  be  obtained  very 
simply  from  the  kd-(3d  equation  obtained  in  3.1.1.  Subsection  3.2  treats  2D  arrays  with  the  dipoles 
parallel  to  the  array  axis.  Since  in  this  case  there  is  no  coupling  between  the  electric  and  magnetic 
dipoles,  this  is  equivalent  to  considering  2D  arrays  of  either  electric  or  magnetic  dipoles,  with  the 
dipoles  parallel  to  the  array  axis. 

In  Section  4  we  obtain  the  kd-(3d  equations  that  determine  the  complex  propagation  constants 
of  traveling  waves  on  3D  arrays  of  lossless  or  lossy  magnetodielectric  spheres.  In  Subsection  4.1 
the  effective  electric  and  magnetic  dipoles  of  the  spheres  are  assumed  to  be  perpendicular  to  each 
other  and  to  the  array  axis,  while  in  Subsection  4.2  the  electric  and  magnetic  dipoles  are  assumed 
to  be  parallel  to  the  array  axis.  Since  in  this  case  there  is  no  coupling  between  the  electric  and 
magnetic  dipoles,  this  is  equivalent  to  considering  3D  arrays  of  either  electric  dipoles  or  magnetic 
dipoles,  with  the  dipoles  parallel  to  the  array  axis. 

In  Section  5  we  obtain  expressions  for  the  effective  constitutive  parameters  of  a  3D  cubic- 
lattice  array  regarded  as  a  homogeneous  medium,  assuming  that  both  \(3d\  -C  1  and  U«1  with 
(3  the  propagation  constant  of  a  transverse  wave  supported  by  the  array.  Two  such  expressions 
are  obtained,  one  from  solving  the  kd  fid  equation,  and  the  other  given  by  the  Clausius-Mossotti 
relations. 

Section  6  deals  with  partially  finite  3D  arrays,  arrays  that  are  finite  in  the  direction  of  the 
array  axis  and  infinite  in  the  transverse  directions.  We  consider  the  field  excited  by  a  plane  wave 
in  the  direction  of  the  array  axis,  incident  from  free  space  on  a  3D  partially  finite  periodic  array 
of  lossless  or  lossy  magnetodielectric  spheres  with  dipoles  perpendicular  to  the  array  axis,  and  on 
a  3D  partially  finite  periodic  array  of  lossless  or  lossy  electric  or  magnetic  dipoles  oriented  parallel 
to  the  array  axis.  Previously  obtained  expressions  for  lossless  dipoles  are  shown  to  be  equally  valid 
for  partially  finite  arrays  of  lossy  dipoles. 

In  Section  7  we  discuss  issues  related  to  the  bidirectionality  and  complex  conjugates  of  waves 
supported  by  periodic  arrays  of  lossless  magnetodielectric  inclusions.  It  is  shown,  consistent  with 
[36],  that  if  (id  is  a  solution  of  a  kd-(3d  equation,  then  so  are  -0d,  ((3d)* ,  and  -(/3 d)*.  However,  for 
ID  and  2D  arrays  it  may  be  necessary  to  choose  alternative  branches  of  the  multivalued  functions 
that  appear  in  the  kd-(3d  equation  if  ((3d)*  and  -((3d)*  are  to  satisfy  the  same  kd-(3d  equation 
satisfied  by  (3d  and  -(3d. 
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Section  8  is  devoted  to  presenting  and  discussing  numerous  kd-fid  diagrams  for  ID,  2D,  and 
3D  periodic  arrays  of  magnetodielectric  spheres  with  equal  permittivity  and  permeability,  diamond 
spheres,  and  silver  nanospheres. 

Appendix  A  contains  an  analysis  of  complex  waves  on  a  2D  slab  of  magnetodielectric  material 
above  a  ground  plane.  Even  though  the  slab  is  a  continuous  structure  and  not  periodic,  considerable 
insight  into  several  of  the  major  issues  of  this  report  can  be  obtained  by  careful  examination  of  the 
slab  problem:  1)  the  derivation  of  the  kd  (id  equations  for  complex  (3  by  analytically  cont  inuing  the 
corresponding  kd  (id  equations  for  real  (3\  2)  the  necessity  of  evaluating  the  multivalued  functions 
that  appear  in  the  kd-(3d  equations  for  ID  and  2D  arrays  using  branches  other  than  the  principal 
one  if  (0d)*  is  to  always  be  a  solution  to  the  same  dispersion  equation  as  (3d  if  the  array  elements 
are  lossless;  and  3)  the  divergence  of  the  complex- wave  fields  obtained  by  direct  summation  of  the 
fields  scattered  from  the  array  dipoles.  Appendix  B  is  devoted  to  a  discussion  of  the  divergence  of 
complex-wave  fields  obtained  by  direct  summation  of  dipole  fields  of  periodic  arrays.  Although  the 
primary  purpose  of  this  report  is  to  derive  and  solve  dispersion  equations  for  waves  in  the  direction  of 
the  array  axis  with  complex  propagation  constants,  such  waves  have  associated  complex  propagation 
constants  in  the  direction  (directions)  transverse  to  the  array  axis.  Accordingly,  in  Appendix  C 
we  present  a  general  discussion  of  complex  waves  on  uniform  and  periodic  waveguides,  and  of  the 
distinguishing  features  of  fast  and  slow  waves. 

In  closing  this  introduction,  we  briefly  discuss  souk*  work  related  to  ours.14  Important  mention 
should  be  made  of  J.  Brown’s  1960  article  [4]  and  R.E.  Collin’s  chapter  on  artificial  dielectrics  [5, 
ch.  12]  consisting  of  periodic  arrays  of  elements  such  as  perfectly  conducting  spheres  and  disks. 
Although  arrived  at  independently  of  Brown  and  Collin,  our  work  can  be  seen  as  a  major  extension 
of  their  work  to  a  full-wave  analysis  of  time-harmonic  dipole  fields  on  periodic  arrays  of  magne¬ 
todielectric  elements.  Although  Brown’s  and  Collin’s  treatment  is  almost  completely  restricted 
to  static  fields,  the  way  they  obtain  the  effective  dielectric  constants  of  artificial  dielectrics  via 
the  interaction  field,  followed  by  an  attempt  to  transform  the  interaction  field  summation  to  a 
rapidly  convergent  form,  parallels  our  method  of  obtaining  the  propagation  constants  of  traveling 
waves  supported  by  periodic  arrays.  Simovski  and  He  [48]  analyzed  a  rectangular  array  of  lossless, 
isotropic  elements  consisting  of  six  0  particles  on  the  faces  of  a  cubic  unit  cell.  Their  procedure 
for  obtaining  the  dispersion  equation  for  the  eigenmodes  of  this  rectangular  lattice  and  for  obtain¬ 
ing  the  effective  parameters  of  a  cubic  lattice  of  the  unit  cells  from  the  solution  to  the  dispersion 
equation,  is  virtually  identical  to  the  procedure  we  developed  independently  in  [20], [21]  and  extend 

14 Despite  its  title  and  its  emphasis  on  the  homogeneous  solution  for  the  propagation  constants  of  waves  on  periodic 
structures,  L.  Brillouin’s  classic  hook  [49]  has  very  little  if  anv  similarity  or  applicability  to  our  work.  It  is  concerned 
for  the  most  part  with  mechanical  rather  than  electromagnetic  interactions  of  particles,  and  with  small  periodic 
perturbations  of  continuous  media.  However,  its  discussion  of  non-rectangular  lattices,  as  well  as  of  propagation 
in  arbitrary  directions  through  periodic  arrays  may  prove  to  be  of  interest  in  possible  extensions  of  our  own  work 
at  a  later  time.  Although  it  might  be  thought  that  the  well-established  science  of  crystallography,  the  science  of 
determining  the  arrangement  of  atoms  in  solids,  with  its  focus  on  periodic  lattices  of  atoms  would  have  an  important 
bearing  on  our  work,  that  is  not  the  case.  The  important  use  of  X-ray  diffraction  techniques  in  crystallography  to 
determine  crystal  structures  is  unrelated  to  the  propagation  of  electromagnetic  waves  through  periodic  arrays  that 
is  the  focus  of  our  work.  In  X-ray  diffraction  each  atom  of  a  crystal  scatters  the  wave  fronts  of  the  incoming  wave 
as  they  arrive  and  so  becomes  surrounded  by  a  series  of  scalar  spherical  wavelets  whose  interference  or  reinforcement 
with  the  spherical  wavelets  from  the  other  atoms  gives  rise  to  a  diffraction  pattern  which  can  then  be  related  to 
the  crystal’s  structure  [50].  Crystal  optics,  the  branch  of  optics  treating  the  behavior  of  light  in  anisotropic  media, 
has  important  applications  to  understanding  and  analyzing  spatial  dispersion:  the  dependence  of  the  electric  and 
magnetic  susceptibility  of  media  such  as  crystals  on  the  propagation  direction  of  light  in  the  media  [51].  Since 
the  analysis  of  3D  arrays  in  this  report  focuses  on  axial  propagation  and  the  conditions  under  which  the  effective 
parameters  derived  from  the  solution  of  the  axial  propagation  dispersion  equations  can  characterize  the  arrays  viewed 
as  isotropic  media,  spatial  dispersion  is  barely  touched  on  here.  However,  spatial  dispersion  is  in  general  an  extremely 
important  topic  in  treatments  of  metamaterials,  and  it  will  be  the  subject  of  forthcoming  work. 
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here  to  lossy  scatterers  and  complex  waves.  The  principal  difference  between  their  work  and  ours 
in  [20], [21]  is  that  they  obtain  the  dispersion  equation  using  an  approximate  method  for  evaluating 
the  infinite  summations  in  the  interaction  field,  whereas  we  use  an  exact,  closed-form  expression  for 
these  summations.  The  work  of  S.  Tretyakov,  especially  [39,  sec.  4.5,  sec.  5.6],  has  been  discussed 
at  some  length  in  [20], [21].  Little  needs  to  be  added  to  that  discussion  here  since  Tretyakov  does 
not  treat  complex  waves  (although  there  is  nothing  in  his  formalism  to  prevent  that),  the  principal 
subject  of  this  report,  or  periodic  arrays  of  magnetodielectric  inclusions.  The  work  of  Alii  and 
Engheta,  in  particular  [43],  has  a  very  strong  family  resemblance  to  our  own  and  has  been  acknowl¬ 
edged  above  in  connection  with  the  analytic  continuation  procedure  to  obtain  complex  waves  that 
plays  a  central  role  in  the  analysis  of  this  report.  The  focus  of  their  work  is  somewhat  different 
from  ours,  however.  Alii  and  Engheta  are  primarily  concerned  with  propagation  along  chains  of 
small  lossy  or  lossless  electrically-polarizable  particles,  not  necessarily  spheres,  rather  than  chains 
of  magnetodielectric  particles,  and  do  not  consider  2D  and  3D  arrays.  Rather  than  emphasize  the 
dependence  of  the  propagation  constant  on  “frequency’1  as  we  do,  they  investigate  how  the  proper¬ 
ties  of  the  propagating  mode  supported  by  the  chain  depend  on  the  polarizability  of  the  particles. 
This  enables  them  to  delineate  regimes  for  leaky-modes  and  guided-modes  that  are  applicable  to 
small  polarizable  particles  of  any  shape,  so  that  their  work  establishes  the  basis  for  using  particle 
shape  as  a  design  variable  if  desired.  Although  not  specifically  concerned  with  traveling  waves  on 
periodic  arrays  of  discrete  scatterers,  A.  Hessel’s  in-depth  mathematically-oriented  treatment  of 
guided  waves  and  traveling  wave  antennas  in  [52,  ch.  10]  is  an  excellent  reference  on  the  general 
subject.  Whereas  our  focus  is  almost  exclusively  on  the  longitudinal  characteristics  of  the  traveling 
waves  that  we  investigate,  Hessel’s  treatment  deals  with  the  properties  of  traveling  waves  in  all 
directions,  transverse  as  well  as  longitudinal,  and  in  so  doing  establishes  the  foundation  for  un¬ 
derstanding  the  waves  in  terms  of  their  complex-plane  representations,  something  that  we  touch 
on  only  briefly  in  Appendices  A  and  C.  The  book  on  traveling  wave  antennas  by  Walter  [53]  has 
some  useful  basic  information  on  fast  and  slow  waves,  and  an  extensive  treatment  of  their  an¬ 
tenna  applications.  Finally,  reference  should  be  made  here  to  photonic  crystals.  Photonic  crystals 
are  optical  artificial  periodic  dielectrics  whose  micro-structure  has  been  engineered  to  control  the 
properties  (propagation,  transmission,  reflection,  band  gaps,  etc.)  of  light  incident  on  the  crystals 
[6]  [8].  Although  the  analysis  of  photonic  crystals  by  casting  Maxwell’s  equations  in  the  form  of 
eigenequations  that  are  solved  computer-numerically  by  sophisticated  discretization  schemes  [7],  [8] 
is  very  different,  from  our  analysis  of  periodic  arrays,  the  generality  of  photonic  crystal  analysis 
can  yield  important  insights  into  matters  of  common  interest  such  as  the  structure  of  dispersion 
diagrams. 

2  ID  ARRAYS 

In  this  section  we  investigate  waves  with  complex  propagation  constants  supported  by  ID  periodic 
arrays  of  lossless  and  lossy  magnetodielectric  spheres.  In  Subsection  2.1  it  is  assumed  that  the 
spheres  can  be  modeled  by  pairs  of  crossed  electric  and  magnetic  dipoles.  In  Subsection  2.2  we 
consider  traveling  waves  supported  by  ID  periodic  arrays  of  lossless  and  lossy  magnetodielectric 
spheres  with  the  dipole  axes  aligned  parallel  with  the  array  axis.  This  reduces  simply  to  traveling 
waves  supported  by  a  ID  array  of  either  electric  or  magnetic  dipoles  aligned  with  the  array  axis 
because  an  electric  (magnetic)  dipole  has  no  radial  or  longitudinal  magnetic  (electric)  field  [38, 
secs.  8.5,  8.6]  and  so  there  is  no  coupling  of  the  electric  dipoles  with  the  magnetic  dipoles  of  arrays 
with  both  the  electric  and  magnetic  dipoles  aligned  with  the  array  axis.  (For  the  same  reason  it 
is  unnecessary  to  consider  ID  arrays  of  electric  and  magnetic  dipoles  with  the  electric  (magnetic) 


9 


dipoles  in  the  direction  of  the  array  axis  and  the  magnetic  (electric)  dipoles  perpendicular  to  the 
array  axis.) 

2.1  ELECTRIC  AND  MAGNETIC  DIPOLES  PERPENDICULAR  TO  THE 
ARRAY  AXIS 

In  this  subsection  we  investigate  the  complex  propagation  constants  of  traveling  waves  supported 
by  ID  periodic  arrays  of  lossless  and  lossy  magnetodielectric  spheres.  It  is  assumed  that  the  spheres 
can  be  modeled  by  pairs  of  crossed  electric  and  magnetic  dipoles  perpendicular  to  each  other,  and 
that  a  spherical  electromagnetic  dipole  field  incident  on  a  sphere  will  excite  an  electric  (magnetic) 
dipole  proportional  to,  and  in  the  direction  of,  the  electric  (magnetic)  field  at  the  sphere  center. 
As  noted  in  the  Introduction,  although  we  refer  in  this  and  the  following  sections  of  the  report  as 
“magnetodielectric  spheres,”  and  assume  that  the  spheres  are  homogeneous,  the  analyses  that  we 
perform  apply  equally  well  to  any  radially  symmetric  array  elements  or  to  elements  that  can  be 
modeled  as  a  pair  of  electric  and  magnetic  dipoles  at  right  angles  to  each  other  and  the  traveling 
wave  direction  provided  that  the  dipoles  are  uncoupled  and  that  an  incident  electric  (magnetic) 
field  at  the  element  center  in  the  direction  of  the  electric  (magnetic)  dipole  excites  only  the  electric 
(magnetic)  dipole  field.  The  radius  of  the  spheres  is  denoted  by  a,  and  the  relative  permittivity  and 
permeability  of  the  spheres  are  denoted  by  er  and  /tr,  respectively,  where  er  and  /ir  are  in  general 
complex.  We  denote  the  separation  between  the  centers  of  adjacent  spheres  by  d,  take  the  2  axis 
to  be  the  axis  of  the  array,  and  assume  an  excitation  of  the  array  with  an  x-directed  electric  field 
and  a  y- directed  magnetic  field.  We  begin  in  exactly  the  same  way  as  we  did  in  our  analysis  of 
traveling  waves  with  real  propagation  constants  supported  by  ID  periodic  arrays  of  lossless  spheres 
[27]  by  obtaining  expressions  for  the  x-directed  electric  field  and  the  y-directed  magnetic  field  at 
the  center  of  the  reference  sphere  at  2  =  0  resulting  from  the  electric  and  magnetic  fields  scattered 
from  all  the  other  spheres  in  the  array.15  In  so  doing  we  make  use  of  the  scattering  equations  [27, 
eq.  31] 


e  h 
90 

CO 

II 

-0 

(2.1a) 

fjO n 

b  -  s  0v 

um,n  urn  y  * 

(2.1b) 

that  relate  the  coefficients  6e,„  (6m,n)  of  the  electric  (magnetic)  dipole  waves  scattered  from  the  nth 
sphere  to  the  incident  electric  (magnetic  field)  E^J  x  (H$/Y0  y)  at  the  center  of  the  nth  sphere. 
In  (2.1)  Se  (Sm)  is  the  normalized  magnetodielectric  sphere  electric  (magnetic)  dipole  scattering 
coefficient.  “Normalized”  means  that  be  n  (bm<n)  is  the  coefficient  of  exp(ifcr )/{kr)  in  the  outgoing 
electric  (magnetic)  dipole  field  in  response  to  the  incident  field  Eq"  x  (Hq"/Yq  y)  at  the  center  of 

IS  All  other  components  cancel  because  of  the  symmetry  of  the  array  lattice. 
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the  x-  ( y -)  directed  electric  (magnetic)  dipole.16  We  then  obtain  the  following  equations 

/ei(/c  +  (3)dn  +  e-i(/3  -  k)dn ' 


W  =  se 


n 


( kdy  + 


i  kd  1  \ 
n  n2 ) 


™  f  i(k  +  (3)dn  _  e— i(/3  -  k)dn 
+  Q (kd)  I  - 


71=  1 


n 


)(m4)’ 


(2.2a) 


and 


(kd)3  =  S„ 


E 

n=l 


si (k  +  f1)dn  _j_  e— i(/3  —  k)dn ' 


n 


1  “  /ei(^  +  fi)dn  _  e-i(/3  -  Ar)dn 

+  ;<“>E  - ; - 


n=l 


where 


<1  = 


0 

bc,o 


(2.2b) 

(2.3) 


In  (2.2)  (5  is  the  real  propagation  constant  of  the  wave  supported  by  the  array. 

In  the  procedure  we  followed  in  our  analysis  of  traveling  waves  with  real  propagation  constants 
supported  by  lossless  spheres  we  were  able  to  obtain  closed-form  expressions  for  the  summations  in 
(2.2).  For  complex  (3  this  is  no  longer  possible  as  in  fact  the  summations  diverge.  Accordingly,  we 
rewrite  (2.2)  in  a  form  that  gives  the  same  result  for  real  (3  as  that  obtained  in  [27]  but  which  can 
be  analytically  continued  into  the  complex  fi  plane.17  To  do  this  we  introduce  the  polylogarithm 
functions  [54],  [55] 


OO  I- 

ZK 


Li AK*)=]CpV’  l2l  <  1 


k=  1 


(2.4) 


which  can  be  analytically  continued  into  the  complex  z  plane  by  means  of  the  iterative  equations 


LiAr(*)  =  j  LljV~lWd t  (2.5a) 

0 

and 

Li^z)  =  -  ln(l  —  z)  .  (2.5b) 

16 As  wo  noted  at  the  beginning  of  this  section  of  the  report,  although  we  refer  to  the  array  elements  as  “magne¬ 
todielectric  spheres,”  our  analysis  applies  equally  well  to  any  array  elements  that  can  be  modeled  as  a  pair  of  crossed 
electric  and  magnetic  dipoles  such  that  an  incident  electric  (magnetic)  field  in  the  direction  of  the  electric  (magnetic) 
dipole  excites  only  the  electric  (magnetic)  dipole  field.  If  the  array  elements  are  homogeneous  spheres  then  Se  and 
Sm  are  the  normalized  Mie  dipole  scattering  coefficients  [27,  eqs.  30a, b], 

Se  —  i-0|  ,  Om  — 

where  fcf  and  of  are  the  electric  and  magnetic  Mie  dipole  scattering  coefficients  defined  in  Stratton  [38,  sec.  9.25, 
eqs.  11,  10].  If  the  array  elements  are  not  homogeneous  magnetodielectric  spheres  then  Sr.  and  Sm  must  be  known 
for  the  results  of  this  and  the  following  sections  of  the  report  to  be  applied. 

17In  Appendix  A  we  justify  this  procedure  by  treating  the  problem  of  determining  the  complex  waves  that  can 
be  supported  by  a  two-dimensional  slab  of  magnetodielectric  material  over  a  perfectly  conducting  ground  plane. 
Although  this  problem  involves  a  continuous  rather  than  a  discrete  periodic  structure,  the  essential  features  of  the 
continuous  and  discrete  periodic  problems  are  identical. 
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Efficient  algorithms  exist  for  calculating  Li/v(z)  for  any  complex  z  (see,  for  example,  [56]).  We  also 
introduce  the  notation 


FN+(kd,  (id)  =  Li*  (e‘(fc  +  W)  +  Liyv  (e_i^  “  W) 

(2.6a) 

FN-(kd,(3d)  =  Li*  (e'(fc  +  W}  -  UN  (e~l(P  ~  fc)d)  . 

(2.6b) 

Then  the  kd  (3d  equations  (2.2)  can  be  written  as 

(kd)3  =  Se  [ 

(kd)2F\+(kd,  (id)  +  i  kdF2+(kd,  (id)  -  Fi+(kd ,  (id) 

+  q  kd(kdFX-(kd,  (id)  +  i  F2-(kd,  /3d))] 

(2.7a) 

(kd)3  =  Sm 

^(kd)2Fi+(kd, (id)  +  i  kdF2+(kd ,  (id )  -  F3+(kd ,  (id) 

+  i  kd(kdF\.(kdJid)  +  i F2_(fcd,/3d))]  . 

(2.7b) 

Letting 

Ei  =  (kd)2F\^.(kd1(3d)  +  i  kdF2+(kd,  (3d)  —  Fz+(kd,  (3d) 

(2.8a) 

and 

E2  =  kd(kdFi-(kd,/id)  +  i  F2_(fcd,/3d)) 

(2.8b) 

in  (2.7)  and  eliminating  q  we 

obtain  the  kd  (id  equation 

(kd)3  -  Se E,  SmE2 

5pE2  (fcd)3  -  5mEi ' 

(2.9) 

As  a  special  case  of  the  kd  (id  equation  (2.9)  for  a  ID  periodic  array  of  lossless  or  lossy  magne¬ 
todielectric  spheres  we  can  obtain  the  kd-/3d  equation  for  a  ID  periodic  array  of  electric  or  magnetic 
dipoles  perpendicular  to  the  array  axis 

(kd)3  =  SE,  (2.10) 

where  S  is  the  normalized  scattering  coefficient  of  either  the  electric  or  the  magnetic  dipoles. 

Using  the  basic  mathematical  relation  that  the  sum  of  logarithms  equals  the  logarithm  of  the 
corresponding  product,  an  alternative  expression  for  Fi+(kd,  (id)  can  be  obtained  from  (2.6a)  with 

N  =  1, 

Fi+(kd,(id)  =  -ikd, —  \n[2(cos  kd  -  cos  (3d)}  .  (2.11) 

This  expression  is  not  equivalent  to  (2.6a)  because  the  branch  cuts  of  (2.6a)  with  N  =  1  in  the 
complex  (id  plane  run  from  kd  to  kd  -I-  ioo  and  from  —kd  to  — kd  —  ioo,  whereas  the  branch  cuts 
of  (2.11)  run  from  kd  to  0  on  top  of  the  real  (id  axis  and  then  to  ioo  along  the  imaginary  (id  axis, 
and  from  -kd  to  0  below  the  real  (id  axis  and  then  to  -ioo  along  the  imaginary  (id  axis  (see  Fig. 
1).  As  a  result,  in  the  regions  1  and  3  of  the  complex  (id  plane  in  Fig.  1,  the  imaginary  part  of 
Fi  +  (kd,(id)  given  by  (2.6a)  is  2n  greater  than  the  value  obtained  with  (2.11). 
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Figure  1:  Branch  cuts  in  the  complex  fid  plane  for  the  two  alternative  expressions  for  F\+(kd,fid). 


2.2  DIPOLES  PARALLEL  TO  THE  ARRAY  AXIS 

In  this  subsection  we  consider  the  complex  propagation  constants  of  traveling  waves  on  ID  periodic 
arrays  of  lossless  and  lossy  magnetodielectric  spheres  with  the  dipoles  parallel  to  the  array  axis.  As 
noted  in  the  introduction  to  this  section,  we  need  only  consider  ID  arrays  of  electric  (or  magnetic 
dipoles)  separately  since  there  is  no  coupling  between  the  electric  and  magnetic  dipoles.  Our 
starting  point  is  the  kd-fid  equation  for  real  fid  [26,  eq.  94] 

OO 

(kd)3  =  2S^ 

n=  1 

where  S  is  the  normalized  dipole  scattering  coefficient.  Instead  of  obtaining  closed-form  expressions 
for  the  summations  as  we  did  in  [26]  we  rewrite  (2.12)  in  terms  of  polylogarithm  functions  so  that 
the  kd-fid  equation  can  be  analytically  continued  into  the  complex  fid  plane 

(kd)3  =  2  S  [-kd.\F2+{kd,  fid)  +  F3+(kd,  fid)}  (2.13) 

where  F2+(kd,0d)  and  F3+(kd,  fid)  are  defined  by  (2.6a),  thereby  obtaining  the  desired  kd-fid 
equation  for  obtaining  the  normalized  complex  propagation  constants  fid. 

3  2D  ARRAYS 

In  this  section  we  consider  traveling  waves  supported  by  2D  periodic  arrays  of  lossless  and  lossy 
magnetodielectric  spheres.  In  Subsection  3.1  it  is  assumed  that  the  spheres  can  be  effectively 
modeled  by  pairs  of  crossed  electric  and  magnetic  dipoles,  each  of  the  dipoles  perpendicular  to 
the  array  axis.  There  are  two  polarizations  of  the  electric  dipoles  to  be  considered,  one  where  the 


ei (k  +  fi)dn  +  e-i (fi  -  k)dn 


(rkdi+i) 


(2.12) 
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dipoles  are  in  the  array  plane  determined  by  the  sphere  centers  and  array  axis,  and  the  other  where 
the  dipoles  are  perpendicular  to  the  array  plane.  These  two  polarizations  are  treated  in  3.1.1  and 
3.1.2,  respectively.  Since  the  electric  and  magnetic  dipoles  of  one  polarization  are,  apart  from  a 
sign  change,  the  magnetic  and  electric  dipoles,  respectively,  of  the  other  polarization,  the  kd-/3d 
equation  for  the  polarization  in  which  the  electric  dipoles  are  perpendicular  to  the  array  plane  can 
be  obtained  very  simply  from  the  kd~0d  equation  for  the  polarization  in  which  the  electric  dipoles 
are  in  the  array  plane.  In  Subsection  3.2  we  consider  traveling  waves  supported  by  2D  periodic 
arrays  of  lossless  and  lossy  magnetodielectric  spheres  with  the  dipole  axes  aligned  parallel  with  the 
array  axis.  This  reduces  simply  to  traveling  waves  supported  by  a  2D  array  of  either  electric  or 
magnetic  dipoles  aligned  with  the  array  axis  because  an  electric  (magnetic)  dipole  has  no  radial  or 
longitudinal  magnetic  (electric)  field  [38,  secs.  8.5,  8.6]  and  so  there  is  no  coupling  of  the  electric 
dipoles  with  the  magnetic  dipoles  of  arrays  with  both  electric  and  magnetic  dipoles  aligned  with  the 
array  axis.  (For  the  same  reason  it  is  unnecessary  to  consider  2D  arrays  of  electric  and  magnetic 
dipoles  with  the  electric  (magnetic)  dipoles  in  the  direction  of  the  array  axis  and  the  magnetic 
(electric)  dipoles  perpendicular  to  the  array  axis.) 


3.1  ELECTRIC  AND  MAGNETIC  DIPOLES  PERPENDICULAR  TO  THE 
ARRAY  AXIS 


3.1.1  ELECTRIC  DIPOLES  IN  THE  ARRAY  PLANE 

We  choose  the  array  axis  to  be  the  2  axis  of  a  Cartesian  coordinate  system  with  equispaced  columns 
of  magnetodielectric  spheres  parallel  to  the  x  axis  located  at  2  =  nd,  n  =  0,  ±1,  ±2,  •  •  • .  In  each 
column  the  spheres  are  centered  at  x  =  rah ,  m  =  0,  ±1,  ±2,  •  *  • .  The  effective  electric  and  magnetic 
dipole  components  of  each  sphere  are  oriented  in  the  x  and  y  direction,  respectively,  so  that  the 
electric  dipoles  lie  in  the  xz  plane,  the  array  plane.  We  assume  an  excitation  of  the  array  with  the 
electric  field  parallel  to  the  x  axis  and  the  magnetic  field  parallel  to  the  y  axis,  and  such  that  all 
the  spheres  in  any  column  of  the  array  are  excited  identically. 

Our  starting  point  is  the  pair  of  equations  [20,  eqs.  8.22a, b]  obtained  by  expressing  the  electric 
and  magnetic  field  at  the  center  of  the  reference  sphere  at  x  =  0,  y  =  0,  2  =  0  as  the  sum  of  the 
electric  and  magnetic  fields  scattered  from  all  the  other  elements  of  the  array, 
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where 
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(3.2) 
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and  Se,  Sm,  and  q  are  as  in  Section  2.  Continuing  just  as  we  did  when  the  normalized  propagation 
constant  f3d  is  real,  from  [20,  eq.  8.25] 
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from  [20,  eq.  8.32] 
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from  [20,  eq.  8.33] 
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and  from  [20,  eq.  8.36] 


l(kh  +  —}=  4  kh  C\2(kh)  +  4  Cl 3{kh)  +  in(kh)2  -  il(kh)3  (3.7) 

i  \  m  J  3 


~  i  khm/  ikh  n 

2  2_. - (  (^l)  - 2  ) 

'  771  V  m  m 

m=l 


=  —  2  ^(fc/i)2ln  2sin^^)  +  khC\2(kh)  +  Cl3(fcJi))  +  i  -  (kh)2  - -(kh)'  (3.8) 

for  0  <  kh  <  2n.  In  (3.4)-(3.8),  //,V 1  is  the  Hankel  function  of  the  first  kind  of  order  n,  Kn  is  the 
modified  Bessel  function  of  order  n,  and  the  Clausen  functions  C\2(kh)  and  Cl :}(kh)  are  defined 
and  given  by  [54] 
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with  the  dilogarithm  and  trilogarithm  functions,  Li2  and  U3,  given  by  (2.4)  and  (2.5). 

To  obtain  analytically  continuable  and  rapidly  convergent  expressions  for  the  Schlomilch  series 
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n=l  n=l 

in  (3.4),  (3.5),  and  (3.6),  following  Linton  [57]  we  introduce  the  notation 
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We  choose  the  branch  cuts  of  \f  [3‘?n  —  1  in  the  complex  /?„,  plane  to  run  from  1  to  1  +  ioo  and  irom 
- 1  to  - 1  -  ioo  with  the  appropriate  branch  of  \J —  1  defined  by 
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Then  from  [57,  eq.  22] 
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where  7,  referred  to  as  C  in  [58],  is  the  Euler  constant  [59,  Table  1.1], 

7  =  0.5772156649- •• 

and 
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with  Bo  and  B2  the  Bernoulli  polynomials  [58,  eq.  9.62] 


B0{x)  =  1,  B-i(x)  =  x2  -  x  +  -  . 


(3.19) 


Truncating  the  series  on  the  RHS  of  (3.14)-(3.16)  at  m  =  10  gives  sufficient  accuracy.  Substituting 
(3.4)-(3.8)  in  (3.1)  and  eliminating  q  from  (3.1a)  and  (3.1b)  we  obtain  the  kd-fid  equation 
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+  4  ^  cos (nfid)  ^ 


Tl=  1 


m=l 


[(27rm)2  +  (kh)2]  K0  (n(d/h)  \/ (2n m)2  -  (kh)2^ 


-  [(27rm)2  -  (kh)2]  K‘2  (n(d/h)  \/ (27rm)2  -  (A:/i)2) 


-  2(kh)2  In 


-  2khC\2(kh)  -  2C\s(kh)  +  i 


(3.23) 


The  rapidly  convergent  analytically  continuable  expressions  (3.14)-(3.16)  are  used  to  calculate 
the  corresponding  Schlomilch  series  in  (3.21)-(3.23).  No  special  treatment  is  necessary  for  the 
Schlomilch  series  involving  the  modified  Bessel  functions  A’o.  K\,  and  A 2  in  (3.21)-(3.23)  since 
even  with  complex  fid  these  series  converge  rapidly  because  of  the  rapid  decay  of  the  modified 
Bessel  functions  with  increasing  n  and  m.  Truncating  these  series  at  n  =  4  and  m  =  4  gives 
sufficient  accuracy. 

As  special  cases  of  the  kd-fid  equation  (3.20)  we  can  obtain  the  kd-0d  equation  for  a  2D  array 
of  electric  dipoles  in  the  plane  of  the  array  and  perpendicular  to  the  array  axis 


(kli)3  =  Se£  1 


(3.24) 


and  the  kd-(3d  equation  for  a  2D  array  of  magnetic  dipoles  perpendicular  to  the  plane  of  the  array 


(kh.)3  =  Sm£3  • 


(3.25) 


3.1.2  ELECTRIC  DIPOLES  PERPENDICULAR  TO  THE  ARRAY  PLANE 

As  shown  in  [20.  sec.  8.2]  the  kd-(id  equation  for  a  2D  array  of  magnetodielectric  spheres  with  the 
electric  dipoles  polarized  perpendicular  to  the  array  plane  can  be  obtained  from  the  kd  (id  equation 
for  a  2D  array  of  magnetodielectric  spheres  with  the  electric  dipoles  in  the  array  plane  simply  by 
interchanging  £1  and  £3.  Thus 
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5^  "  (kh)3  -  5m£  1  ] 

where  £j,  £2,  and  £3  are  given  (3.21)-(3.23). 

As  special  cases  of  the  kd  (id  equation  (3.26)  we  can  obtain  the  kd-0d  equation  for  a  2D  array 
of  magnetic  dipoles  in  the  plane  of  the  array  and  perpendicular  to  the  array  axis 

(kh)3  =  Sm£i  (3.27) 

and  the  kd-(3d  equation  for  a  2D  array  of  electric  dipoles  perpendicular  to  the  plane  of  the  array 

(kh)3  =  5e£3  .  (3-28) 

3.2  DIPOLES  PARALLEL  TO  THE  ARRAY  AXIS 

As  we  noted  in  the  introduction  to  this  section,  there  is  no  need  to  consider  arrays  with  both  electric 
and  magnetic  dipoles  parallel  to  the  array  axis  since  there  is  no  coupling  between  the  electric  and 
magnetic  dipoles.  Without  loss  of  generality  we  can  consider  2D  arrays  of  electric  dipoles  parallel 
to  the  array  axis.  Our  starting  point  is  the  kd  (id  equation  [20,  eq.  7.21]  obtained  by  expressing 
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the  electric  field  incident  on  the  reference  electric  dipole  at  the  origin  as  the  sum  of  electric  dipole 
waves  scattered  from  all  the  other  electric  dipoles  in  the  array 
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for  0  <  kh  <  2n  with  the  Clausen  functions  CI2 {kh)  and  Cl3(fc/i)  given  by  (3.9).  Substituting 
(3.31)  and  (3.32)  in  (3.29)  we  obtain  the  computational  form  of  the  kd-fid  equation 
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+  [(27rm)2  -  (kh)2]  K 2  (n(d/h)  y/^hrw^p-  -  (kh)2) 


-2( 

The  Schlomilch  series  with  Hq^  and  H^]  are  evaluated  using  the  analytically  continuable  expres¬ 
sions  (3.14)  and  (3.16).  The  Schlomilch  series  with  A'o  and  K2  can  be  truncated  at  n  =  4  and 
m  =  4  with  sufficient  accuracy. 
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4  3D  ARRAYS 

In  this  section  we  investigate  waves  with  complex  propagation  constants  supported  by  3D  periodic 
arrays  of  lossless  and  lossy  magnetodielectric  spheres.  In  Subsection  4.1  it  is  assumed  that  the 
spheres  can  be  effectively  modeled  by  pairs  of  crossed  electric  and  magnetic  dipoles,  each  of  the 
dipoles  perpendicular  to  the  array  axis.  In  Subsection  4.2  we  consider  traveling  waves  supported 
by  3D  periodic  arrays  of  lossless  and  lossy  magnetodielectric  spheres  with  the  dipole  axes  aligned 
parallel  with  the  array  axis.  This  reduces  simply  to  traveling  waves  supported  by  a  3D  array  of 
either  electric  or  magnetic  dipoles  aligned  with  the  array  axis  because  an  electric  (magnetic)  dipole 
has  no  radial  or  longitudinal  magnetic  (electric)  field  [38,  secs.  8.5,  8.6]  and  so  there  is  no  coupling 
of  the  electric  dipoles  with  the  magnetic  dipoles  of  arrays  with  both  the  electric  and  magnetic 
dipoles  aligned  with  the  array  axis.  (For  the  same  reason  it  is  unnecessary  to  consider  3D  arrays 
of  electric  and  magnetic  dipoles  with  the  electric  (magnetic)  dipoles  in  the  direction  of  the  array 
axis  and  the  magnetic  (electric)  dipoles  perpendicular  to  the  array  axis.) 

4.1  ELECTRIC  AND  MAGNETIC  DIPOLES  PERPENDICULAR  TO  THE 
ARRAY  AXIS 

In  this  subsection  we  consider  traveling  waves  supported  by  3D  periodic  arrays  of  lossless  and 
lossy  magnetodielectric  spheres.  It  is  assumed  that  the  spheres  can  be  modeler!  by  pairs  of  crossed 
electric  and  magnetic  dipoles,  each  of  the  dipoles  perpendicular  to  the  array  axis.  The  analysis 
performed  here  is  equally  applicable  to  any  3D  periodic  arrays  whose  elements  can  be  modeled 
by  a  pair  of  crossed  electric  and  magnetic  dipoles  at  right  angles  to  each  other  such  that  only  an 
incident  electric  (magnetic)  field  at  the  element  center  in  the  direction  of  the  electric  (magnetic) 
dipole  excites  only  the  electric  (magnetic)  dipole  field.  We  choose  the  array  axis  to  be  the  z 
axis  of  a  Cartesian  coordinate  system  with  equispaced  planes  of  magnetodielectric  spheres  normal 
to  the  z  axis  located  at  z  =  nd.,n  =  0,  ±1,±2,  ••••  In  each  plane  the  spheres  are  centered  at 
x  —  mh,y  =  lh,l,m  =  0,  ±1,  ±2,  •  •  • .  The  effective  electric  and  magnetic  dipole  components  of 
each  sphere  are  oriented  in  the  x  and  y  direction,  respectively.  We  assume  an  excitation  of  the 
array  with  the  electric  field  parallel  to  the  x  axis  and  the  magnetic  field  parallel  to  the  y  axis,  and 
such  that  all  the  spheres  in  any  column  of  the  array  are  excited  identically. 

Our  starting  point  is  the  pair  of  equations  [20,  9.28], [21,  eqs.  35a,b]  obtained  by  expressing  the 
electric  and  magnetic  field  incident  on  the  reference  sphere  at  x  =  0,  y  =  0,  z  =  0  as  the  sum  of  the 
electric  and  magnetic  dipole  waves  scattered  from  all  the  other  elements  of  the  array 
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and  Se,  Sm,  and  q  are  as  in  Section  2.  Proceeding  just  as  we  did  in  our  analysis  of  traveling  waves 
with  real  propagation  constants  supported  by  3D  arrays  of  lossless  magnetodielectric  spheres,  from 
[20,  eq.  9.30], [21,  eq.  37] 
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oo  oo 
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(4.4) 


with  the  Clausen  functions  Cl2  ami  CI3  given  by  equations  (3.9),  and  with  0  <  kh  <  27r.  From  [20, 
eq.  9.76], (21,  eq.  41] 
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Substituting  (4.4)  and  (4.5)  in  (4.1),  using  [58,  eq.  8.521(1)] 
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and  eliminating  g  from  (4.1a)  and  (4.1b)  we  obtain  the  kd-(3d  equation 
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and 


E2  =  27r  kh 


sin  [3d 


cos  /id  —  cos  fcd 


-  4tt  kh  E  sin(n/W)  E  e~n(d/'l)  +  /2)  -  (*h)2 


(4.9) 


71=1 


m,/=  -  to 
(m.l)*<0.0> 


In  (4.6)-(4.9),  Jo  and  Vo  are  the  ordinary  Bessel  functions  of  the  first  and  second  kind.  This  kd-/3d 
equation  is  analytically  continuable  into  the  complex  [M  plane  as  it  stands  simply  by  using  the 
cosine  and  sine  of  a  complex  argument.  Truncating  the  series  in  Ei  and  E2  with  the  decaying 
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exponential  and  with  A'0  at  n  =  m  =  l  =  4  yields  sufficient  accuracy.  The  series  T.YQ(lkh)  in  £1  is 
evaluated  by  [58,  eq.  8.522(3)] 


E«“)~  i(7  +  l"£)-2£ 

/=1  ^  '  1=1 

with  7  the  Euler  constant  given  by  (3.17).  Truncating  the  series  on  the  RHS  of  (4.10)  at  l  =  10 
gives  sufficient  accuracy  for  our  purposes. 

As  special  cases  of  the  kd-0d  equation  (4.7)  we  can  obtain  the  kd-(3d  equation  for  a  3D  array 
of  electric  or  magnetic  dipoles  perpendicular  to  the  array  axis 

(kli)3  =  51! i  .  (4.11) 

where  5  is  the  normalized  scattering  coefficient  of  either  the  electric  or  the  magnetic  dipoles. 


^(2 ttI)2  -  (Jfch)2  27rZ 


,  0  <  kh  <  27r 


(4.10) 


4.2  DIPOLES  PARALLEL  TO  THE  ARRAY  AXIS 

In  this  subsection  we  consider  the  complex  propagation  constants  of  traveling  waves  on  3D  periodic 
arrays  of  lossless  and  lossy  magnetodielectric  spheres  with  the  dipoles  parallel  to  the  array  axis.  As 
noted  in  the  introduction  to  this  section,  we  need  only  consider  3D  arrays  of  electric  (or  magnetic 
dipoles)  separately  since  there  is  no  coupling  between  the  electric  and  magnetic  dipoles.  Our 
starting  point  is  the  kd-/3d  equation  [20,  eq.  7.21] 


(kh.f  =  5  ei 


i  nfid 


n~  —  oc 

n*  0 


£ 

m./=— oo 


gi  khpmin 
Pmln 


—  2i 

Pmln 


kh  + 


i\  (nd/h 

Pmln )  Pmln 


+ 


i kh  1  \ 

(kh)2  +  - -j—  ) 

\  Pmln  Pmln ' 


TTl 


2  +  Z2 


~  eikhpmio 

+  /  - 

“  PmlO 

m,i=-oc 

(m.l)^(O.O) 


Pmln  /  P min 


( ,  .o  i  kh  1 

(  (^0  +  “ - 

\  PmlO  PmlO  > 


(4.12) 


where  5  is  the  normalized  electric  dipole  scattering  coefficient, 

Pmln  =  \J  m2  +  l2  +  (nd/h)2 

and 
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From  (20,  eq.  7.64] 
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(4.14) 
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and  from  [20,  eq.  7.71] 
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Substituting  (4.14)  and  (4.15)  in  (4.12)  and  using  (4.6)  and  [20,  eqs.  B.13] 
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(4.15) 

(4.16) 

(4.17) 
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In  (4.15)-(4.18),  J„,  Yn,  Kn,  and  //r(,1),  are  the  nth  order  ordinary  Bessel  functions  of  the  first  and 
second  kind,  the  modified  Bessel  function,  and  the  Hankel  function  of  the  first  kind,  respectively, 
and  CI2  and  CI3  are  the  Clausen  functions  given  by  (3.9).  The  kd-ftd  equation  is  analytically 
continuable  into  the  complex  (3d  plane  as  it  stands  simply  by  using  the  cosine  of  a  complex  argument. 
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Truncating  the  series  in  Q  with  the  decaying  exponential  and  with  A'o  at  n  =  m  =  l  =  4  yields 
snfficient,  accuracy.  The  series  Sin (/ kh )  is  evaluated  by  (4.10)  while  t  he  series  —  F2 (Ikh )  is  evaluated 
by  [20,  eqs.  B.14  -  B.15],[21,  eq.  A9] 
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and  Bq  and  B2  are  the  Bernoulli  polynomials  given  by  (3.19).  Truncating  the  series  on  the  RHS  of 
(4.19)  at  m  =  10  gives  sufficient  accuracy. 


5  EFFECTIVE  PERMITTIVITY  AND  PERMEABILITY 

When  the  magnetodielectric  spheres  of  a  3D  cubic-lattice  array  are  sufficiently  close  so  that  both 
ftd.  <C  1  and  kd  <  1  where  (3  is  the  real  propagation  constant  of  a  transverse  traveling  wave 
supported  by  the  array,  we  argue  in  [20,  sec.  9.3], [21,  sec.  3.2]  that  the  array  can  be  regarded 
macroscopically  as  a  homogeneous  medium,  which  we  refer  to  as  the  array  medium.  (See,  however, 
Footnote  23  below.)  The  array  medium  is  characterized  by  an  effective  or  bulk  relative  permittivity 
and  effective  relative  permeability  nf  that  determine  the  propagation  constant  of  a  traveling 
wave  in  the  direction  of  the  array  axis  perpendicular  to  the  orientations  of  the  crossed  electric 
and  magnetic  dipoles  by  which  the  spheres  are  modeled.18  We  then  show  how  and  nf  can  be 
obtained  from  the  parameters  available  to  us  in  solving  the  kd-(3d  equation.19  For  complex  (3d  the 
assumption  that.  (3d  <  1  is  replaced  by  the  assumption  that  \(3d\  <  1,  the  final  expressions  for  the 
effective  constitutive  parameters  are  the  same  as  for  real  (3d ,  and  the  derivation  of  these  expressions 
differs  only  slightly  from  that  of  the  derivation  for  real  (3d.  Accordingly  we  can  omit  most  of  the 
details  here. 

18In  general  the  array  medium  is  anisotropic  and  and  //', ''  do  not  determine  the  propagation  of  waves  traveling 
in  directions  other  than  along  the  array  axis.  If  the  array  elements  are  radially  symmetric  magnetodielectric  spheres 
then  the  directions  of  the  effective  electric  and  magnetic  dipoles  of  the  array  are  established  by  the  traveling  wave, 
and  as  the  number  of  spheres  per  unit  volume  becomes  large  (kd  <?C  1  as  well  as  fid  1)  the  array  medium  becomes 
increasingly  isotropic.  If,  however,  the  directions  of  the  effective  electric  and  magnetic  dipoles  of  the  array  elements 
are  fixed  independently  of  the  traveling  wave,  as  they  are  for  split-ring  resonators  for  example,  then  the  array  medium 
is  anisotropic  no  matter  how  closely  spaced  the  elements. 

19 At  the  time  of  writing  [20]  and  [21]  we  were  unfortunately  unaware  that  in  their  2003  paper  [48],  Simovski 
and  He  had  developed  a  method  for  obtaining  the  effective  parameters  from  the  dispersion  equation  for  a  traveling 
wave  supported  by  a  cubic-lattice  array  of  magnetodielectric  elements  (cubic  unit  cells  with  an  fi-shaped  perfectly 
conducting  particle  on  each  face)  that  is  essentially  the  same  as  ours.  The  principal  difference  between  their  method 
and  ours  is  that  their  expressions  for  the  effective  parameters  are  based  on  a  form  of  the  dispersion  equation  obtained 
by  approximating  the  infinite  summations  involved,  whereas  our  method  is  based  on  the  exact  dispersion  equation 
obtained  using  closed-form  expressions  for  the  infinite  summations. 
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As  in  [20,  eq.  9. 93], [21,  eq.  56]  the  starting  point  is  the  relation,  valid  for  transverse  plane 

20 

waves 

yA  =  v/Wff  «?ff  (5-1) 

where,  as  throughout,  we  shall  assume  that  the  sign  of  the  square  root  is  chosen  to  keep  the  real 
part  of  (3  positive  in  the  exp(i (3z)  dependence.  [20,  eq.  9. 100], [21,  eq.  62]  is  replaced  by 


m  _  pf  -  1 

p  —  1  /i‘r?fT 


where  c  is  the  speed  of  light,  and 

m 

—  =cq 
V 


as  before,  where  q  is  defined  by  (2.3)  and  obtained  from  (4.7). 
From  (5.3)  and  (5.2)  we  obtain 


jeff  _  1  ^eff  9 


(5.2) 


(5.3) 


(5.4) 


Equations  (5.1)  and  (5.4)  form  a  pair  of  simultaneous  equations  which  can  be  solved  for  the  two 
unknowns  and  Letting  1Z  =  fid/{kd)  we  obtain  as  in  [20,  eq.  9.107,  9. 108], [21,  eq.  68] 

eft  _  +  l)  eff  _  ^(1  +^9)  Ac  c \ 

*r  "  1  +  lZq  ’  ^  ~  7 ZTq  ‘  (5‘5) 

It  should  be  noted  that  if  the  opposite  choice  of  the  sign  of  the  square  root  in  (5.1)  is  made,  then 
(3d,  7 7,  and  q  change  sign  so  that  e*ff  and  /z*ff  given  by  (5.5)  remain  the  same.21 

In  [20,  sec.  9. 2], [21,  sec.  3.2]  we  also  obtain  expressions  for  the  bulk  permittivity  and  perme¬ 
ability  of  the  array  medium  based  on  the  Clausius-Mossotti  relation  [35,  sec.  8-1],  [60,  sec.  2-4], 
assuming  that  (3d  1  and  kd  C  1,  and  that  the  transverse  spacing  of  the  array  elements  h  equals 

the  spacing  d  in  the  direction  of  the  array  axis.  These  expressions  are  derived  with  no  reference  to 

20For  a  3D  array  of  electric  or  magnetic  dipoles  parallel  to  the  array  axis  (direction  of  propagation),  Maxwell’s 
equations  show  that  the  array  is  not  characterized  by  a  homogeneous,  isotropic,  scalar  permittivity  and  permeability 
satisfying  5.1  and  5.2. 

21  It  might  appear  from  the  expressions  (5.5)  for  the  bulk  parameters  that  for  lossless  arrays,  as  the  light  line  is 
approached  by  the  dispersion  curve  so  that  77  %  1  (and  f and  ^  become  reciprocals),  both  and  /x?ff  should 
approach  the  value  of  +1.  However,  a  little  thought  and  numerical  examples  such  as  [20,  figs.  37  and  38;  figs. 
39  and  40], [21,  figs.  25  and  26;  figs.  27  and  28]  show  that  when  the  light  line  is  approached  by  a  backward- wave 
(negative-slope)  segment  of  the  kd-pd.  diagram,  both  c?ff  and  are  negative.  Furthermore,  numerical  results  show 
that  and  ^  approach  the  value  of  -1  as  the  light  line  is  approached  by  a  backward- wave  segment  of  the  dispersion 
curve,  only  for  array  spheres  of  equal  relative  permittivity  and  permeability.  For  array  spheres  of  unequal  relative 
permittivity  and  permeability,  approaches  a  value  different  from  -1  (and,  of  course,  ficTn  approaches  the  reciprocal 
value).  How  is  it  that  the  computer- program  implementation  of  (5.5)  gives  results  that  are  so  much  at  variance  with 
what  one  might  expect  from  a  cursory  examination  of  (5.5)?  The  answer  is  simply  that  77  and  q  are  not  independent 
of  one  another  and  if  the  fine  numerical  details  of  their  behavior  are  examined  as  the  light  line  is  approached  by  a 
back  ward- wave  segment  of  the  kd-Pd  diagram,  the  apparent  contradiction  disappears.  If,  for  example,  we  consider 
the  case  of  [20,  figs.  37  and  38], [21,  figs.  25  and  26]  where  the  relative  permittivity  and  permeability  of  the  array 
spheres  are  equal,  then  it  is  found  from  computer  results  that  as  the  light  line  is  approached  by  the  backward-wave 
segment,  q  very  nearly  approaches  the  value  of  -1  and  77  %  1  +  c,  |c|  1  so  that  77(77  +  q)  zz  e,  l+77q%  — c, 

and  hence  «  -1.  A  similar  analysis  can  be  performed  for  the  case  of  array  spheres  with  unequal  relative 

permittivity  and  permeability. 


the  solution  of  the  kd-/3d  equation  and  are  equally  valid  here  provided  that  \fid\  1  and  kd  l:22 


.off, CM 
cr 


2D  +  3 
3  -  B 


(5.6) 


when1 


eff,CM 

Mr 


2  A  +  3 

3  -  ,4 


J3  =  - 


6nitf 

(kd)* 


67ria]c 

(kd)* 


(5.7) 

(5.8) 

(5.9) 


and  6®c  and  a*]c  are  the  Mie  electric  and  magnetic  dipole  scattering  coefficients,  respectively,  given 
in  [38,  sec.  9.25,  eqs.  11,10]. 

The  derivations  given  in  this  section  assume  that  \fid]  and  H  «  1.  In  practice,  however,  as 
mentioned  in  the  Introduction,  these  conditions  can  be  relaxed  somewhat  to  the  “rule  of  thumb” 
that  \/3d\  and  kd  be  less  than  about  1  [21,  sec.  3.2]. 23 


6  PARTIALLY  FINITE  ARRAYS 

In  [20,  sec.  11], [21,  sec.  5]  we  obtain  expressions  for  the  field  at  any  point  on  the  array  axis  (other 
than  at  elements  of  the  array)  excited  by  a  plane  wave  incident  from  free  space  on  a  3D  partially 
finite  periodic  array  of  lossless  magnetodielectric  spheres.  The  array  is  finite  in  the  direction  of  the 

22The  conditions  needed  in  general  to  derive  the  C lausi us- M ossot t i  effective  parameter  expressions  are  that  there 
should  be  many  dipoles  in  a  macroscopic  volume  in  which  the  phase  and  magnitude  of  the  dipoles  is  approximately 
uniform.  In  the  context  of  this  report,  these  conditions  take  the  form  \fid\  1  and  kd  1  which  in  practice  can  be 
relaxed  to  \f3d\  <  1  and  kd  <  1. 

23Simovski  [61]  has  shown  that  even  if  \0d\  and  kd  are  both  «:  1  it  is  still  possible  that  the  effective  parameters 
obtained  from  solving  the  dispersion  equation  for  the  propagation  constant  of  a  transverse  traveling  wave  in  the  axial 
direction  of  a  cubic-lattice  array  may  not  accurately  characterize  the  array  regarded  as  an  isotropic,  homogeneous 
medium.  If,  as  for  example  in  the  first  backward  branch  of  the  dispersion  diagram  shown  in  [22,  fig.  2],  \@d\  is 
changing  very  rapidly  within  a  very  small  kd  interval  (frequency  band),  then  even  a  small  change  in  the  dispersion 
diagram  that  can  be  expected  if  propagation  in  a  direction  other  than  axial  is  considered,  can  result  in  a  large 
change  in  the  associated  values  of  the  effective  permittivity  and  permeability.  Examining  [22,  figs.  3,4],  for  example, 
where  the  effective  parameters  corresponding  to  the  dispersion  diagram  of  [22,  fig.  2]  are  plotted  versus  kd,  it  is 
seen  that  the  effective  parameters  are  changing  extremely  rapidly  with  kd  in  the  vicinity  of  the  kd  interval  when 
they  are  negative.  Thus,  the  frequency  interval  where  \0d\  is  small  and  the  effective  parameters  are  both  negative, 
may  no  longer  be  associated  with  negative  effective  parameters  if  propagation  in  a  non-axial  direction  is  considered. 
It  should  be  emphasized,  however,  that  the  extremely  narrow  frequency  band  of  the  dispersion  diagram  backward 
branch  in  the  example  of  [22]  we  have  referred  to  here,  is  by  no  means  typical,  and  for  wider  frequency  bands  with 
slower  variation  of  \0d\  with  frequency,  the  effective  parameters  obtained  from  the  solution  to  the  axial  propagation 
dispersion  equation  can  be  expected  to  reasonably  well  characterize  the  array  regarded  as  an  isotropic,  homogeneous 
medium  when  both  \0d\  and  kd  are  small  .  One  further  important  point  should  be  noted  here.  The  Clausius-Mossotti 
expressions  for  the  effective  parameters  are  derived  with  no  reference  to  the  propagation  direction.  Consequently, 
the  frequency  regions  where  the  effective  parameters  obtained  from  the  solution  of  the  axial-propagation  dispersion 
equation  can  lx?  expected  to  be  a  good  characterization  of  the  array  regarded  as  an  isotropic  medium,  are  such  that 
these  effective  parameters  agree  reasonably  well  with  the  Clausius-Mossotti  effective  parameters.  Thus  a  combined 
plot  of  the  dispersion  equation  effective  parameters  and  the  Clausius-Mossotti  effective  parameters  is  a  very  useful 
tool  for  estimating  the  frequency  regions  where  the  array  can  be  expected  to  behave  as  a  quasi-isotropic  medium, 
and  conversely  where  isotropy  cannot  be  expected  to  hold.  In  the  example  of  [22]  cited  in  this  footnote,  it  is  seen 
that  the  plots  of  the  dispersion  equation  effective  parameters  increasingly  diverge  from  the  Clausius-Mossotti  effective 
parameter  plots  as  the  doubly-negative  kd  interval  is  approached,  a  clear  indication  that  isotropy  may  not  hold  in 
the  doubly-negative  kd  region. 
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array  axis  and  infinite  in  the  directions  transverse  to  the  array  axis.  The  direction  of  incidence  of 
the  illuminating  plane  wave  is  parallel  to  the  array  axis,  normal  to  the  interface  between  the  array 
and  free  space.  The  spheres  are  modeled  by  pairs  of  crossed  electric  and  magnetic  dipoles,  each  of 
the  dipoles  perpendicular  to  the  array  axis.  The  expressions  obtained  for  the  field  on  the  array  axis 
in  [20,  sec.  11], [21,  sec.  5]  are  completely  independent  of  the  kd-(3d  equation  and  are  equally  valid 
as  they  stand  for  a  3D  partially  finite  periodic  array  of  lossless  or  lossy  magnetodielectric  spheres. 

The  corresponding  expression  for  the  field  on  the  array  axis  of  a  3D  partially  finite  periodic 
array  of  electric  dipoles  parallel  to  the  array  axis,  applicable  to  either  lossless  or  lossy  dipoles,  is 
given  in  [20,  sec.  11.2].  The  same  expression  is  valid  as  well  for  a  3D  partially  finite  array  of 
magnetic  dipoles  parallel  to  the  array  axis. 


7  BIDIRECTIONALITY  AND  COMPLEX  CONJUGATES  OF 


WAVES 


In  this  section  we  discuss  the  bidirectionality  and  complex  conjugates  of  the  propagation  constants 
of  traveling  waves  in  the  direction  of  the  array  axis  supported  by  ID,  2D,  and  .‘ID  periodic  arrays 
of  lossless  scatterers.  Our  starting  point  is  the  paper  [36]  (also  see  Appendix  C)  in  which  it  is 
proved  that  a  reciprocal  (lossy  or  lossless)  waveguide  (uniform  or  periodic)  that  supports  a  traveling 
wave  with  propagation  constant  ft  also  supports  a  corresponding  traveling  wave  with  propagation 
constant  —ft,  and  that  for  every  traveling  wave  with  a  complex  propagation  constant  ft  on  a  lossless 
reciprocal  uniform  or  periodic  waveguide  there  exists  a  corresponding  traveling  wave  with  complex 
conjugate  propagation  constant  ft* .  Regarding  bidirectionality,  it  is  clear  from  the  analysis  we 
have  presented  in  Sections  2,  3,  and  4,  either  from  inspection  of  the  starting  summations  or  of 
the  final  forms  of  the  dispersion  equations,  that  if  ftd  is  a  solution  of  the  dispersion  equation  for  a 
ID,  2D,  or  3D  array,  then  so  is  -ftd.  More  complicated,  however,  is  the  relation  of  our  results  to 
the  above  complex  conjugate  propagation  constant  theorem.  For  the  lossless  3D  arrays  considered 
in  this  report,  it  is  easy  to  show  that  if  ftd  is  a  complex  solution  to  the  dispersion  equation  then 
so  is  (ftd)*.  For  the  lossless  ID  and  2D  arrays  we  have  considered,  however,  numerical  results 
show  that  (ftd)*  is  not  necessarily  a  solution  to  the  same  dispersion  equation  satisfied  by  complex 
ftd.  However,  in  these  ID  and  2D  cases,  (ftd)*  is  a  solution  to  the  dispersion  equation  when  the 
multivalued  functions  that  appear  in  the  dispersion  equation  are  evaluated  with  branches  other 
than  the  principal  one.  We  first  look  at  the  3D  arrays,  and  then  discuss  the  2D  and  ID  arrays. 

Consider  the  kd-ftd  equation  (4.7)  for  the  3D  array  of  magnetodielectric  spheres  with  the 
dipoles  oriented  perpendicular  to  the  array  axis.  Dividing  the  LHS  and  RHS  of  (4.7)  by  Se  and 
Sm ,  respectively,  we  obtain 


*  E2  (kh)3 


(7.1) 


Referring  to  (4.8)  let 


V,  =E,+i?(fc/i)3. 


(7.2) 


Then,  using  the  relations  [20,  eq.  4. 61], [21,  eq.  45] 


(7.3) 
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where  i/>e  and  ipm  are  the  phases  of  the  normalized  lossless  Mie  scattering  coefficients  Se  and  Sm, 
respectively,  we  see  that 

(kh)3  ^  2/ilx3  cost/>< 


and 


Se 

(kh)3 


-s t  =  =(*fc)3  -e; 

3  sint/>e 


E  2 

3V  '  sin  4’m  1 


(7.4a) 

(7.4b) 


so  that  the  dispersion  equation  (4.7)  is  equivalent  to 

i(MOa 

^  Sin^e  
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i(fc/»)3 
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sin  ipn 
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Referring  to  (4.8)  and  (4.9)  we  see  that  changing  (3d  to  {/3d)*  changes  both  sides  of  (7.5)  to  their 
complex  conjugates,  so  that  if  (3d  is  a  solution  of  (7.5)  then  so  is  {(3d)*.  Exactly  the  same  line  of 
reasoning  applied  to  a  3D  array  of  dipoles  parallel  to  the  array  axis  shows  that  if  (3d  is  a  solution 
of  the  kd (3d  equation  (4.17)  then  so  is  {/3d)*. 

Next,  let  us  consider  2D  arrays.  Just  as  investigating  the  propagation  of  waves  on  a  2D  slab  in 
Appendix  A  gives  much  insight  into  the  analytic  continuation  procedure  that  plays  a  central  role 
in  the  main  body  of  the  report,  so  our  treatment  of  the  slab  can  also  give  us  considerable  insight 
into  a  wave  with  propagation  constant  {(3d)*  corresponding  to  one  with  propagation  constant  (3d. 
The  kd-(3d  equation  derived  in  Appendix  A  for  a  TE  wave  supported  by  a  2D  slab  is 


tan  *yd  —  7 


7  d 


Hblod 


with 


and 


')d=[{kd)2-m2\x'2 

7o d  =  [(k0d)2  -  (Pd)2] 1/2 
(kd)2  —  u/2/lf  =  (ko)Vrfr  • 


(7.6) 

(7.7a) 

(7.7b) 

(7.7c) 


If  ( (3d )*  is  substituted  for  Qd  in  (7.6),  tailed),  7 d,  and  7 0d  all  change  into  their  complex  conjugates, 
but  the  factor  of  i  in  the  denominator  of  the  RHS  of  (7.6)  remains  unchanged.  Hence  (fid)*  is  not 
a  solution  of  (7.6)  as  it  stands.  However,  if  the  alternative  branch  for  the  square  root  function 
defining  70 d  is  used,  then  instead  of  70 d  in  (7.6)  we  have  —  7o<i  and  (fid)*  is  then  a  solution  of  the 
modified  form  of  (7.6).  Much  the  same  idea  applies  to  the  kd-fid  equation  (3.20)  for  a  2D  array  of 
magnetodielectric  spheres  with  the  electric  dipoles  in  the  plane  of  the  array.  The  term  i  \/ /3q  —  1 
appears  in  the  expressions  (3.14),  (3.15)  and  (3.16)  for  calculating  the  Schlomilch  series  that  appear 
in  Ei,  E2,  and  E3  given  by  (3.21)-(3.23).  Referring  to  (3.11) 


v  7o 


(7.8) 


(note  that  k,  the  free-space  wave  number  in  the  main  body  of  the  report  is  identical  to  k0  in 
Appendix  A).  Exactly  analogous  to  the  way  (fid.)*  is  a  solution  of  the  slab  dispersion  equation 
(7.6)  if  the  alternative  branch  for  the  square  root  function  defining  70  is  used,  so  if  -y/fio  —  1  is 
substituted  for  \Jfil  -  1  in  (3.14),  (3.15),  and  (3.16),  and  these  expressions  then  substituted  in 
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(3.20),  then  (/3d)*  is  found  to  be  a  solution  of  the  modified  form  of  (3.20)  when  /3d  is  a  solution 
of  the  original  form  of  (3.20).  The  same  conclusion  holds  for  the  kd-(3d  equation  (3.26)  of  the  2D 
array  of  magnetodielectric  spheres  with  the  electric  dipoles  perpendicular  to  the  plane  of  the  array, 
and  the  kd  (3d  equation  (3.33)  of  a  2D  array  of  electric  dipoles  parallel  to  the  array  axis.  Although 
in  general  this  report  does  not  consider  propagation  in  the  direction  transverse  to  the  array  axis, 
for  the  2D  arrays  considered  in  this  report  propagation  in  the  direction  y  normal  to  the  plane  of 

the  array  is  seen  to  be  given  by  e^0^  =  e^  ”  @  This  wave  is  proper  or  improper  if  3(70) 

is  positive  or  negative,  respectively. 

In  the  case  of  ID  arrays,  consideration  of  waves  supported  by  a  magnetodielectric  rod  does  not 
give  any  insight  into  why  (3 *  is  often  not  a  solution  of  the  same  dispersion  equation  satisfied  by 
/3  unless  alternate  branches  of  the  dispersion  function  are  taken  into  account,  unlike  the  situation 
for  the  2D  slab  that  we  have  just  discussed.  (For  the  magnetodielectric  rod,  it  is  easy  to  see  from 
the  expressions  for  the  dispersion  equations  for  TE,  TM,  HE,  and  HM  modes  given  in  [62,  sec. 
3.2]  that  if  (3  is  a  solution  of  the  dispersion  equation  then  so  is  /3*.24)  However,  just  as  in  the  2D 
continuous  and  periodic  cases,  if  (3d  is  a  solution  of  the  kd-(3d  equation  for  the  ID  array,  then  so  is 
(/3d)*  provided  that  branches  of  the  multivalued  functions  other  than  the  principal  ones  are  used 
when  the  kd  (3d  function  is  evaluated. 

Refer  to  the  kd-(3d  equation,  (2.9),  for  a  ID  array  of  magnetodielectric  spheres  with  dipoles 
oriented  perpendicular  to  the  array  axis.  In  this  equation  £i  and  £2  are  defined  in  terms  of  sums 
and  differences  of  the  polylogarithm  functions  Li jv(z),  N  =  1, 2, 3.  To  obtain  higher  order  branches 
of  the  multivalued  functions  that  appear  in  £1  and  £2  we  need  expressions  for  the  higher  order 
branches  of  Liyv(^)  given  the  principal  values  of  these  functions.  Since  Lii(z)  =  —  ln(l  —  2),  higher 
order  branches  of  Li  1(2)  are  obtained  simply  by  adding  integer  multiples  of  27ri  to  the  principal 
value  of  -ln(l  —  2).  To  obtain  the  relation  of  the  higher  order  branches  of  Li2 (2)  and  Li3(2)  to 
their  principal  values  we  proceed  as  follows.  The  poly  logarithm  functions  satisfy  the  differential 
recursion  relation  [63,  eq.  18] 

dLi.s(2)  / *7  -% r\\ 

Z— ^  =  Lis_,(z).  (7.10) 

Let 

Li„(z)  =  Li„(z)prjncipaj  +  2mmf(n)  ln"-1(z)  .  (7.11) 

Then  from  (7.10)  f(n)  satisfies  the  recursion  relation 

/(n  +  1)  =  >  /(l)  =  1-  (7-12) 

Then  /( 2)  =  1,  /( 3)  =  and 

Li2(z)  =  Li2(z)prjncipai  +  2rmriln(z),  m  =  0,  ±1,  ±2,  •  •  •  (7.13a) 

and 

Li3(z)  =  Li3(z)prindpai  +  mm  ln2(z),  m  =  0,  ±1,  ±2,  •  •  •  .  (7.13b) 

24  As  an  example,  the  dispersion  equation  for  TE  modes  [62,  eq.  3.2-1 7a]  is 

J\  (ha)  _  h  i  {qa) 

ha  Jo)  ha)  qaKo(qa) 

where  a  is  the  rod  radius,  Jo  and  J\  are  the  ordinary  Bessel  functions  of  order  0  and  1,  respectively,  Kp  and  K\  are 
the  modified  Bessel  functions  of  order  0  and  1,  respectively,  ha  =  y/n2(k0a)2  -  (0a)2  ,qa  =  y/(n2  -  1  )(kpd)2  -  (ha)2, 
and  n  =  y/Crfir  is  the  refractive  index  of  the  rod  material.  Changing  3  to  3*  has  the  effect  of  changing  both  sides  of 
(7.9)  to  their  complex  conjugates  so  that  if  (7.9)  is  satisfied  by  ft  it  is  also  satisfied  by  ft*. 


Since  the  logarithm  function  in  (7.13)  itself  has  branches,  (7.13)  can  be  more  completely  written 
as 


Li2(z)  =  Li2 ( -)principal  +  2m7ri  ,n(2)  ““  4m/c7r2,  m  =  0,  ±1,  ±2,  •  •  •  ,  k  =  0,  del,  ±2,  •  •  •  (7.14a) 


and 


Li:*(z)  =  Li.3(2)principaj  +  m7ri  In 2(z)  -  imkn2  ln(2)  -  4rafc27r3i,  m  =  0,  ±1,  ±2,  •  •  •  , 

k  =  0,  ±1,  ±2,  •  •  •  .  (7.14b) 

Equation  (7.14a)  is  identical  to  [55,  eq.  3.13].  Using  (7.14)  it  is  then  a  simple  matter  to  evaluate 
(2.9)  for  any  choice  of  branches.  In  practice  the  branches  given  by  (7.13)  have  been  found  sufficient 
to  show  numerically  that  if  fid  is  a  solution  of  (2.9)  evaluated  with  the  principal  branches  of 
the  multivalued  functions  then  (fid)*  is  also  a  solution  of  (2.9)  provided  branches  other  than  the 
principal  can  be  chosen.  Of  course  the  same  conclusion  holds  for  the  kd-(3d  equation,  (2.13),  for  a 
ID  array  of  dipoles  parallel  to  the  array  axis. 

In  Appendix  C  it  is  shown  that,  for  lossless  reciprocal  waveguides,  the  only  difference  between 
the  full  propagation  behavior  (that  is,  combined  longitudinal  and  transverse  propagation)  of  the 
(3  and  —fi*  waves  decaying  in  the  positive  z  direction  is  the  sign  of  the  phase.  This  result  follows 
from  the  general  theorem,  derived  from  Maxwell’s  equations,  which  states  that  for  every  solution 
(E,  H)  on  a  lossless  reciprocal  waveguide  there  is  another  solution  (E*,  -H*)  [36].  Also,  since  it 
can  be  shown  [36],  as  discussed  at  the  beginning  of  this  section,  that  if  a  reciprocal  waveguide 
has  a  complex-wave  solution  with  propagation  constant  in  the  2  direction  (3 ,  then  there  exists  a 
corresponding  solution  with  propagation  constant  in  the  z  direction  —0.  Therefore,  it  is  of  interest 
to  know  when  and  why  ( 0d )*  is  a  solution  of  the  same  principal- branch  ID  and  2D  kd  fid  equation 
as  is  satisfied  by  0d ,  and  when  it  is  not.  As  to  “when?”,  a  considerable  number  of  numerical 
examples  support  (but  do  not  necessarily  establish)  the  conclusion  that  in  the  fast-wave  region 
(3d)*  is  never  a  solution  of  the  principal-branch  dispersion  equation  satisfied  by  0d.  (It  is,  as 
stated  above,  always  a  solution  of  the  dispersion  equation  evaluated  with  alternative  branches  of 
the  multivalued  functions  involved.)  In  the  slow- wave  region,  (0d)*  is  generally  a  solution  of  exactly 
the  same  principal-branch  dispersion  equation  satisfied  by  fid,  apart  from  some  examples  where 
the  magnitude  of  the  imaginary  part  of  a  solution  fid  is  very  large,  or  when  5R (fid)  is  very  close  to 
the  light  line.  We  are  unable  at  the  present  time  to  support  these  conclusions  analytically.  One 
possible  conjecture  based  on  the  physics  is  that  in  the  fast-wave  region,  solutions  of  the  kd-fid 
equations  for  our  ID  and  2D  arrays  of  magnetodielectric  spheres  correspond  to  improper  “leaky 
waves”  that  are  not  part  of  the  complete  set  of  modes  and  thus  are  not  all  found  as  solutions  to 
the  principal-branch  dispersion  equation.  The  slow-wave  solutions,  on  the  other  hand,  generally 
correspond  to  proper  waves  that  are  part  of  the  complete  set  of  modes  and  thus  it  can  be  expected 
that  all  possible  solutions  are  determined  by  the  principal-branch  dispersion  equation,  except  where 
3(/3d)  is  too  large  or  'Sl(fid)  is  very  close  to  kd  (the  light  line). 

Tamir  and  Oliner  [64]  (see  Appendix  C)  have  looked  into  the  question  of  whether  physically 
realistic  sources  can  produce  a  significant  amount  of  energy  in  the  four  complex  waves  decaying 
in  the  positive  2  direction,  the  two  proper  complex  waves  (outgoing  and  incoming)  and  the  two 
improper  complex  waves  (outgoing  and  incoming).  They  find  that  the  two  proper  complex  waves 
can  be  appreciably  excited  by  realistic  sources.  In  addition,  in  the  near- field  region  of  waveguides, 
the  outgoing  improper  complex  wave  (the  classic  leaky  wave)  can  dominate  the  fields.  However,  the 
incoming  improper  complex  wave  cannot  apparently  be  excited  by  finite-magnitude  sources  of  finite 
extent.  In  particular,  if  fi  refers  to  a  proper  complex  wave,  then  both  the  fi  and  -fi*  waves  can 
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be  excited  to  an  appreciable  extent  by  finite  sources.  However,  if  (3  refers  to  an  improper  outgoing 
wave  (leaky  wave),  it  can  be  appreciably  excited  by  realistic  sources  but  the  corresponding  —  /3* 
improper  incoming  wave  cannot. 

8  NUMERICAL  RESULTS 

In  this  section  we  present  a  variety  of  kd-0d  diagrams  obtained  with  high-accuracy  computer 
implementations  of  the  expressions  derived  in  the  previous  sections.  Our  objective  in  this  section 
is  primarily  to  display  a  representative  selection  of  kd  -fid  diagrams  without  attempting  to  give 
physical  explanations  of  the  behavior  shown  in  these  diagrams.  Alu  and  Engheta  in  a  series  of  papers 
[65], [66], [43], [67]  have  made  an  important  contribution  towards  obtaining  physical  understanding  of 
the  behavior  of  dispersion  diagrams  of  periodic  arrays  of  electric  dipoles  in  terms  of  nanocircuit  and 
nanotransmission  lines.  Further  work  along  these  lines  is  needed  for  arrays  of  magnetodielectric 
spheres  for  which  both  electric  and  magnetic  dipoles  and  their  interactions  must  be  taken  into 
account. 

The  examples  of  dispersion  diagrams  we  show  here  are  for  the  most  part  designed  to  demon¬ 
strate  what  happens  when  the  dispersion  diagrams  obtained  in  our  previous  papers  and  reports 
[27], [20], [21]  under  the  assumption  that  /3d  is  real,  are  extended  to  allow  for  complex  f3d.  Some 
important  notations  and  conventions  used  in  our  plots  of  the  dispersion  diagrams  are  as  follows.  We 
consistently  use  solid  and  dashed  lines  to  plot  3t(0d)  and  0(/?d),  respectively,  and  use  dotted  lines 
to  plot  light  lines  (kd  =  $t(0d)).  Harmonic  time  dependence  of  exp(-icjf),  u  >  0  is  assumed.  In 
the  description  of  the  arrays,  a  denotes  the  radius  of  the  magnetodielectric  spheres,  and  d  denotes 
the  distance  between  the  centers  of  adjacent  array  elements.  We  plot  3i(0d)  in  the  interval  [0, 7r] . 
Thus  the  sign  of  <3(0d)  corresponding  to  this  positive  (+z  direction)  phase  velocity  determines 
whether  the  wave  is  a  “forward”  [9(/3d)  >  0]  or  “backward”  [^(0d)  <  0]  complex  wave  [17],  If  we 
required  that  ^s(/3d)  always  be  positive,  then  3 i(0d)  would  be  negative  for  “backward”  waves  and 
%t(0d)  would  have  to  be  plotted  in  the  interval  [ — 7r,  7 r]  instead  of  in  the  interval  [0,  tt]. 

The  group  velocity  [cd(kd)  /  d3R(0d)\  is  not  necessarily  in  the  direction  of  decay  of  the  complex 
wave  because  the  proof  of  the  equality  between  group  and  energy-transport  velocities  requires  that 
0  be  real  [68].  Moreover,  for  proper  complex  waves  on  lossless  arrays,  the  total  energy  transport 
across  an  xy  plane  (z  =  constant  plane)  is  zero.  Also,  the  proof  of  the  equality  between  group  and 
energy-transport  velocities  applies  to  continuous  media  and  thus  requires  not  only  that  0  be  real, 
but  also  that  0d  -C  1  and  H<1.  Therefore,  in  reference  to  complex  waves,  the  terms  “forward” 
and  “backward”  simply  denote  that  the  phase  velocity  and  exponential  decay  are  in  the  same  or 
opposite  directions,  respectively,  and  not  necessarily  that  the  phase  and  group  velocities  are  in  the 
same  and  opposite  directions.  For  the  portions  of  the  kd-  0d  diagrams  where  0  is  real,  the  wave 
corresponding  to  a  positive  (negative)  slope  of  the  kd  3d  diagram  is  often  referred  to  as  a  “forward 
(“backward”)  wave.  However,  in  this  paper  we  will  reserve  the  terms  “forward”  and  “backward” 
wave  solely  for  complex  waves  as  defined  above,  and  use  positive  and  negative  group  velocity  as 
synonymous  with  positive  and  negative  slope  of  the  portions  of  the  k  0  diagram  where  0  is  real. 

As  we  have  noted,  the  polylogarithm  functions  in  the  ID  array  kd-0d  equations,  and  the  square 
root  functions  in  the  2D  array  kd-0d  equations  are  multivalued.  Obviously,  time  and  space  do  not 
allow  us  to  calculate  and  show  the  dispersion  diagrams  for  all  choices  of  the  branches  of  these 
multivalued  functions.  Accordingly  we  show  kd-0d  diagrams  corresponding  to  a  choice  of  branches 
other  than  the  principal  branches  of  these  multivalued  functions  for  only  a  small  selection  of  arrays. 
(See  especially  the  kd-0d  diagrams  in  Subsubsections  8.3.1  and  8.3.2  for  ID  and  2D  arrays  of  silver 
nanospheres  with  the  electric  dipoles  oriented  parallel  to  the  array  axis.  Also  see  Mode  3  in  Fig. 
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4.)  As  noted  in  the  Introduction,  because  of  the  difficulty  of  finding  complex  zeroes  of  the  kd-(id 
equation  in  a  given  region  of  the  complex  (3d  plane,  it  is  possible  for  all  the  arrays  considered  that 
there  are  curves  additional  to  those  shown  in  the  kd (id  diagrams  of  this  section  of  the  paper. 

8.1  MAGNETODIELECTRIC  SPHERE  ARRAYS 

8.1.1  LINEAR  (ID)  ARRAYS  OF  MAGNETODIELECTRIC  SPHERES 

In  this  subsubsection  we  consider  twro  linear  arrays,  one  with  er  =  f.ir  —  10,  and  the  other  with 
€r  =  /^r  =  20.  For  a  linear  array  of  magnetodielectric  spheres  with  er  =  /xr  =  10,  a/d  =  0.45,  and 
the  dipoles  transverse  to  the  array  axis,  the  kd-/3d  diagram  for  (id  real  is  given  in  [27,  fig.  12].  In 
Fig.  2  we  show  the  first  branch  of  this  diagram  with  a  positive-slope  segment  for  kd  from  0.724 
to  0.884  in  which  (3d  increases  from  the  light  line  to  n  at  kd  =  0.884  ( ka  =  0.398)  close  to  the 
resonance  frequency  of  kd  =  0.900  (ka  =  0.405),  followed  by  a  negative-slope  segment  for  kd  from 
0.884  to  0.959  where  (3d  decreases  from  n  to  the  light  line.  The  extension  of  this  diagram  into  the 
fast-wave  region  is  shown  in  Fig.  3.  This  extension  has  two  distinct  components.  In  the  upper  one 
(kd  >  0.928),  the  backward- wave  5R(/3d)  portion  of  Fig.  2  crosses  the  light  line,  decreases  to  0  at 
kd  =  0.960,  and  then  increases  again  as  a  forward  wave  to  the  light  line  at  kd  =  0.998  where  it 
ends.  In  this  portion  of  the  kd  (3d  diagram  Q((3d)  decreases  from  0  to  —1.28  as  9 1((3d )  decreases 
from  the  light  line  to  0  (the  backward-wave  segment),  then  changes  sign  and  increases  to  2.26  as 
3?(/?d)  increases  to  the  light  line  (the  forward-wave  segment).  The  end  of  the  plots  of  3?(/?d)  and 
3( fid )  at  kd  =  0.998  deserves  comment.  This  is  not  the  end  of  the  entire  kd-(3d  diagram.  As  kd 
continues  to  increase  a  new  branch  of  the  kd  ( id  diagram  begins  at  kd  %  1.22  (we  do  not  show  this) 
but  in  the  kd  interval  from  0.998  to  1.22  we  have  been  unable  to  find  any  discrete  wave  at  all.  This 
stop-band  interval  beginning  at  the  light  line  is  seen  repeatedly  in  other  ID  and  2D  array  figures 
of  this  section  of  the  paper.  Unless  our  numerical  search  routine  for  complex  zeroes  of  the  kd-(3d 
equation  has  overlooked  some  solutions,  it  appears  that  in  these  frequency  ranges  of  ID  and  2D 
figures  with  no  solutions  of  the  kd-(3d  equation,  the  only  fields  that  can  exist  are  the  “space  waves” 
from  the  continuous  spectra  corresponding  to  integrals  in  the  complex  (3d  plane  along  the  branch 
cuts  of  the  kd  (3d  equation  (see  Appendix  A).  These  space  waves  require  external  sources  for  their 
excitation.  If,  for  example,  a  ID  array  of  magnetodielectric  spheres  is  excited  by  an  external  source, 
then  in  addition  to  discrete  traveling  waves  that  can  channel  some  of  the  incident  energy  along  the 
array,  there  will  also  be  a  continuous  spectrum  of  scattered  fields  in  other  directions,  the  space 
waves.  In  certain  frequency  ranges  external  sources  can  excite  only  these  space  waves.  These  are 
exactly  the  frequency  ranges  where  there  are  no  solutions  to  the  kd  (3d  equation.  (For  3D  arrays  the 
kd  (3d  equations  do  not  contain  functions  with  branch  cuts  so  that,  there  are  only  discrete  modes. 
Consequently  none  of  the  3D  kd  (3d  diagrams  display  stop  bands  beginning  at  the  light  line.)  The 
lower  component  of  the  extension  of  the  kd  (id  diagram  of  Fig.  2  into  the  fast-wave  region  begins 
with  a  backward- wave  segment,  at  first  near  vertical  until  it  reaches  the  3?(/3d)  curve  of  Fig.  2 
with  which  it  is  then  almost  identical  for  a  narrow  frequency  range.  The  negative-slope  segment 
then  takes  off  from  the  5?(/3d)  curve  of  Fig.  2  at  kd  &  0.872,  and  decreases  to  0  at  kd  =  0.916, 
following  which  there  is  a  forward-wave  segment  that  ends  at  the  light  line  at  kd  «  0.930.  What 
is  particularly  interesting  here  is  that  there  is  a  small  section  of  the  backward-wave  segment  in 
the  slow-wave  region  where  (3d  is  complex,  a  transition  region  referred  to  in  [69]  as  “the  winding 
down  of  the  leaky-wave  solution.”  25  Despite  the  similarity  between  this  backward-wave  segment  of 

25The  structure  analyzed  in  [69]  is  very  different,  from  ours,  being  a  dielectric-filled  rectangular  waveguide  with  an 
asymmetric  slit  in  its  top  wall  and  with  an  air-filled  parallel  plate  stub  guide  above,  and  the  frequency  behavior  is 
the  opposite  of  ours. 
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Fig.  3  and  the  transition  region  of  [69,  figs.  3,4]  we  do  not  want  to  assert  that  the  backward-wave 
segment  in  the  fast -wave  region  of  Fig.  3  corresponds  to  a  leaky  wave;  that  is,  an  outward-energy- 
propagating  improper  complex  wave.2*’  In  this  paper  we  consider  only  the  longitudinal  propagation 
constants  of  waves  propagating  in  the  direction  of  the  array  axis  and  do  not  attempt  in  general  to 
investigate  the  behavior  of  these  waves  in  the  direction(s)  transverse  to  the  array  axis.  In  Fig.  4  we 
show  the  kd  (id  diagram  for  the  same  array  but  with  additional  modes  displayed,  demonstrating 
the  complexity  of  the  kd  (id  diagram  for  complex  waves.  Mode  3  in  Fig.  4  was  obtained  using 
(2.11)  instead  of  (2.6a)  for  calculating  F\  +  (kd,  (Id,). 

As  an  example  to  demonstrate  the  importance  that  both  magnetic  and  electric  dipoles  play  in 
the  structure  of  the  dispersion  diagram,  we  show  the  kd-/3d  diagrams  corresponding  to  Figs.  2 
and  3  when  magnetic  dipoles  are  omitted  in  the  calculations.  The  Mie  electric  dipole  scattering 
coefficient  was  computed  as  before  but  the  Mie  magnetic  dipole  scattering  coefficient  was  set  equal 
to  zero.  In  Fig.  5  we  show  the  dispersion  diagram  obtained  when  (id  is  restricted  to  be  real  and  in 
Fig.  6  we  show  the  dispersion  diagram  obtained  when  0d  is  allowed  to  be  complex.  Comparison  of 
corresponding  figures  shows  the  great  difference  between  the  diagrams  obtained  with  and  without 
the  presence  of  magnetic  dipoles.  The  (id  curve  of  Fig.  5  terminates  at  the  value  of  (id  =  n  with 
a  negative-slope  segment  from  (id  =  1.203  (kd  =  0.900)  to  (id  =  n  (kd  =  0.884).  In  Fig.  6  we  see 
that  for  kd  <  0.900  the  dispersion  diagram  is  double-valued,  the  initial  positive-slope  section  with 
0!(/ id )  =  0,  and  a  section  with  3?(/?d)  =  7 r  and  S(/3d)  decreasing  from  2.49  at  kd  =  0.800  to  0  at 
kd  =  0.884.  The  upper  part,  kd  =  0.935  to  0.972,  of  the  extension  of  the  dispersion  diagram  into  the 
fast-wave  region  with  5R(/?d)  increasing  from  0  to  the  light  line  is  very  similar  to  the  corresponding 
segment  in  Fig.  3.  As  in  Fig.  3  there  is  no  continuation  of  the  kd-(id  diagram  after  $t((id)  reaches 
the  light  line.  The  lower  part  of  the  extension  into  the  fast-wave  region  is  a  backward-wave  segment 
that  begins  in  the  slow-wave  region  at  5R(/?cf)  =  1.230  (kd  =  0.900)  starting  at  a  point  on  the  lower 
branch  of  the  dispersion  curve  and  then  crosses  the  light  line  into  the  fast-wave  region.  As  in  Fig. 
3  this  results  in  there  being  a  small  interval  where  $i((id)  is  in  the  slow-wave  region  and  3(/ id )  is 
not  zero,  a  behavior,  as  noted  in  our  discussion  of  Fig.  3,  similar  to  the  transition  region  of  [69, 
figs.  3,4]  with  the  small  portion  of  this  segment  in  the  slow-wave  region  referred  to  as  “the  winding 
down  of  the  leaky-wave  solution.’' 

For  a  linear  array  of  magnetodielectric  spheres  with  er  =  fir  =  20,  a/d  =  0.45,  and  the  dipoles 
transverse  to  the  array  axis,  the  first  branch  of  the  kd-(id  diagram  for  (id  real  is  given  in  [27, 
fig.  15]  which  we  reproduce  here  as  Fig.  7.  The  extension  of  this  diagram  for  complex  (id  is 
shown  in  Fig.  8.  The  principal  features  of  Fig.  8  are  very  similar  to  those  of  Fig.  3  for  the  ID 
er  =  fir  =  10  array,  allowing  for  the  lower  frequency  region  of  Fig.  8.  The  lower  branch  of  the 
diagram  goes  to  n  at  ka  =  0.211  (kd  =  0.469)  close  to  the  resonance  frequency  of  the  spheres  at 
ka  =  0.214  (kd  =  0.475).  The  most  interesting  feature  is  the  portion  of  the  diagram  for  kd  between 
0.476  and  0.481,  the  backward-wave  segment  for  1 R((id)  that  takes  off  from  the  first  branch  of  the 
diagram,  crosses  the  light  line,  and  continues  into  the  fast-wave  region  with  ?R((id)  decreasing  to 
0.  There  is  again  a  small  section  of  this  segment  in  the  slow- wave  region  where  (id  is  complex,  the 
transition  region  referred  to  in  [69]  as  “the  winding  down  of  the  leaky-wave  solution.”  In  Fig.  9  we 
show  the  kd-(id  diagram  for  the  same  array  but  with  additional  modes  displayed.  The  essential 
features  of  this  figure  are  similar  to  those  of  Fig.  4.  As  in  Fig.  4,  Mode  3  in  Fig.  9  was  obtained 
using  (2.11)  instead  of  (2.6a)  for  calculating  F\+(kd,  (id). 

For  a  linear  array  of  magnetodielectric  spheres  with  er  =  /ir  =  20,  a/d  —  0.45,  and  the  electric 
dipoles  parallel  to  the  array  axis,  the  kd  (id  diagram  for  fid  real  is  given  in  [20,  fig.  41]  which  we 
reproduce  here  as  Fig.  10.  The  extension  of  this  diagram  into  the  fast-wave  region  is  shown  in  Fig. 

26 A  detailed  treatment  of  leaky  waves  and  leaky- wave  antennas  can  be  found  in  [37]. 
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11.  For  kd  between  0.407  and  0.453,  fid  is  complex  with  dt(fid)  describing  a  fish  hook  from  the 
light  line  into  the  fast- wave  region  and  then  back  again  to  the  light  line.  Note  that  the  bottom  of 
the  fish  hook  is  not  a  continuation  of  a  slow-wave  region  curve  into  the  fast-wave  region  but  begins 
at  the  light  line.  Note  also  that  this  portion  of  the  fish  hook  is  a  negative-slope  wave  but  3( fid, )  is 
positive.  For  kd  between  0.453  and  0.484,  fid  is  real,  increasing  from  the  light  line  to  7 r  close  to  the 
resonance  frequency  of  the  spheres  at  ka  =  0.214  (kd  =  0.475).  For  values  of  kd  >  0.484,  there  is  a 
strongly  attenuated  wave  with  5 R(fid)  remaining  equal  to  n  and  Q(/?d)  increasing  rapidly. 

8.1.2  2D  ARRAYS  OF  MAGNETODIELECTRIC  SPHERES 

For  fid  real,  the  kd-fid  diagram  for  a  2D  array  of  magnetodielectric  spheres  with  e.r  =  fir  =  20, 
a/d  =  0.45,  and  the  electric  dipoles  parallel  or  perpendicular  to  the  plane  of  the  array,  is  given  in 
[20,  fig.  36], [21,  fig.  24]  and  is  reproduced  here  as  Fig.  12.  Its  extension  for  complex  fid  is  shown 
in  Fig.  13.  The  lower  branch  of  the  diagram  goes  to  n  at  kd  =  0.449  (ka  =  0.202)  fairly  close  to 
the  resonance  frequency  of  the  spheres  at  kd  =  0.475  (ka  =  0.213)  while  the  negative-slope  branch 
of  the  diagram  starts  at  kd  =  0.472  (ka  =  0.212).  The  bandgap  between  kd  =  0.4487  and  0.4721  in 
Fig.  12  is  seen  to  correspond  to  an  attenuated  wave  with  IR(fid)  =  7 r  and  3( fid )  increasing  from  0.0 
to  1.35  at  the  center  of  the  bandgap.  The  behavior  of  the  dispersion  diagram  from  kd  =  0.486  to 
0.495  is  especially  interesting,  and  we  show  a  detail  of  this  region  in  Fig.  14.  We  have  encountered 
behavior  like  this  before  in  Figs.  3,  6,  and  8,  where  a  fast-wave  region  upper  segment  with  'Si(fid) 
increasing  from  0  to  the  light  line  is  joined  to  the  curve  of  !R(/3d)  in  the  slow-wave  region  by  a 
backward-wave  segment  which  includes  a  small  transition  interval  in  the  slow-wave  region  where 
fid  is  complex.  Note  that,  as  in  Figs.  3,  6,  and  8,  there  is  a  cut-off  of  the  kd-fid  diagram  at  the 
upper  frequency  of  the  plot.  In  Fig.  15  we  show  the  kd-fid  diagram  for  the  same  array  but  with  an 
additional  mode  displayed.  Note  the  high  attenuation  of  this  mode  (Mode  2)  compared  with  that 
of  Mode  1.  Figure  16  is  a  more  detailed  version  of  Fig.  15. 

8.1.3  3D  ARRAYS  OF  MAGNETODIELECTRIC  SPHERES 

For  fid  real,  the  kd-fid  diagram  for  a  3D  array  of  magnetodielectric  spheres  with  er  =  [ir  =  20, 
a/d  =  0.45,  and  the  electric  and  magnetic  dipoles  normal  to  the  array  axis  is  given  in  [20,  fig. 
37], [21,  fig.  25]  and  is  reproduced  here  as  Fig.  17.  When  complex  solutions  of  the  kd  fid  equation 
are  allowed  we  obtain  the  kd-fid  diagram  of  Fig.  18.  Mode  1  is  identical  to  the  dispersion  curve 
of  Fig.  17,  and  an  additional  mode  with  high  attenuation  is  shown,  also  obtained  as  a  solution  to 
the  kd-fid  equation.  Mode  1  goes  to  ir  at  kd  =  0.450  (ka  =  0.203)  close  to  the  resonance  frequency 
of  the  spheres  at  kd  =  0.475  (ka  =  0.213).  This  figure  demonstrates  that  complex  modes  can  be 
supported  by  lossless  3D  structures.  No  total  power  is  carried  by  these  modes,  however. 

In  Fig.  19  we  show  the  dispersion  diagrams  for  the  principal  mode  of  Fig.  18  when  the 
permeability  of  the  spheres  has  a  loss  tangent  of  0.002  and  0.02.  The  real  part  of  the  dispersion 
diagrams  are  so  close  to  the  dispersion  diagram  when  the  loss  tangent  is  zero  that  even  though 
the  three  curves  are  plotted  in  Fig.  19  it  is  impossible  to  distinguish  them.  The  permeability  loss 
tangents  of  0.002  and  0.02  do,  however,  result  in  an  imaginary  part  of  the  propagation  constant 
fid  as  shown  in  the  dashed  line  for  the  0.002  loss  tangent  and  the  dot-dashed  line  for  the  0.02  loss 
tangent.  For  the  0.002  loss  tangent,  the  imaginary  part  of  fid  results  in  a  power  loss  per  wavelength 
[  -201og10(exp(-27r|3(/3d)|/A:d))  =  -20(27r\^(fid)\/kd)  lnlO  ]  of  between  2.04  and  4.96  dB  per 
wavelength  in  the  negative-slope  segment  of  the  dispersion  diagram,  while  for  the  0.02  loss  tangent, 
the  imaginary  part  of  fid  results  in  a  proportionally  greater  maximum  power  loss  of  between  20.1 
and  50.1  dB  per  wavelength.  For  a  linear  array  with  er  =  5ft (fir)  —  10  (not  shown),  the  power  loss 


for  a  0.002  permeability  loss  tangent,  lies  between  1.00  and  1.34  clB  per  wavelength,  and  for  a  0.02 
permeability  loss  tangent  the  power  loss  lies  between  10.0  and  13.4  dB  per  wavelength.  Currently 
we  are  not  aware  of  any  magnetodielectric  material  at  frequencies  above  a  few  hundred  MHz  with 
both  relative  permittivity  and  permeability  having  real  parts  reasonably  close  to  the  same  value 
of  about  10  or  greater,  and  a  permeability  loss  tangent  sufficiently  small,  that  can  be  used  for  the 
fabrication  of  low-loss,  isotropic  doubly  negative  (DNG)  materials.  (For  values  of  the  real  parts  of 
the  permittivity  and  permeability  appreciably  less  than  10,  it  does  not  appear  possible  to  obtain 
solutions  of  the  kd-(3d  equation  with  both  \(3d\  and  kd  <  1  as  required  for  the  array  to  reasonably 
approximate  a  DNG  homogeneous  isotropic  medium.) 

8.1.4  EFFECTIVE  PARAMETERS  FOR  3D  ARRAYS  OF  MAGNETODIELEC¬ 
TRIC  SPHERES 

In  Figs.  20  and  21  we  show  the  kd  [3d  diagrams  and  plots  of  the  real  and  imaginary  parts  of  the 
corresponding  effective  parameters,  and  for  two  3D  cubic-lattice  arrays  of  magnetodielectric 
spheres:  1)  a  lossless  array  with  er  =  (13.8,0.0)  and  pr  =  (11.0,0.0);  and  2)  a  lossy  array  with 
er  =  (13.8,0.1)  and  //r  =  (11.0,0.0).  In  both  arrays  the  dipoles  are  perpendicular  to  the  array  axis 
and  a/d  —  0.4.  (3d  is  complex  for  the  lossless  as  well  as  for  the  lossy  array.  For  clarity  the  plots  in 
both  figures  have  been  limited  to  the  range  0.7  <  kd  <  0.9.  In  the  real  (3d  dispersion  diagram  for 
the  first,  array  there  are  stop-bands  in  the  intervals  kd  =  [0.7893,0.8078]  and  kd  =  [0.8738,0.8812] 
where  no  wave  with  a  real  (3  can  propagate.  In  Fig.  20,  in  these  frequency  intervals  (3d  is  complex  so 
that  the  wave  attenuates.  In  the  interval  kd  =  [0.7893,0.8078]  %t((3d)  =  n  and  ^f((3d)  increases  from 
0  to  0.301  and  then  decreases  to  0.  In  the  interval  kd  =  [0.8738,0.8812],  5J(/3d)  =  0  and  |3f(/3d)| 
increases  from  0  to  0.091  and  then  decreases  to  0  and  thus  is  too  small  to  be  seen  in  the  scale  of  the 
plot  (as  are  the  corresponding  imaginary  parts  of  the  effective  parameters).  In  the  first,  and  second 
positive-slope  branches  of  the  kd  (3d  diagram,  kd  =  [0.7, 0.7893]  and  [0.8738, 0.9],  it  is  seen  from  Fig. 
21  that  e®1*  and  /z®ff  are  real  and  positive,  while  in  the  negative-slope  branch,  kd  =  [0.8078,0.8738], 
the  effective  parameters  are  both  negative,  an  example  of  how  a  DNG  medium  can  be  formed  from 
a  periodic  array  of  magnetodielectric  spheres.  Although  in  the  interval  kd  =  [0.7893,0.8078],  \(3d\ 
is  too  large  for  the  array  to  behave  as  a  continuous,  homogeneous,  isotropic  medium,  nevertheless 
viewed  formally  the  complex  values  of  the  effective  parameters  are  consistent  with  the  complex 
value  of  (3d  in  this  interval  and  can  yield  a  reflection  coefficient  for  a  slab  of  the  array  that  is  a 
reasonable  approximation  to  the  actual  reflection  coefficient  .  A  consequence  of  the  formal  extension 
of  the  bulk  parameters  beyond  the  region  where  they  are  physically  meaningful  is,  however,  that 
the  imaginary  parts  of  dr’lf  and  can  become  negative  in  this  interval. 

The  effect  of  loss  on  the  behavior  of  K(e<;fr)  is  seen  from  Fig.  20  to  be  significant  only  in  the 
bandgap  region  kd  —  [0.7893,0.8078].  Whereas  tor  the  lossless  array,  5?(e'r^)  has  a  positive  and 
negative  infinite  asymptote  at,  kd  =  0.80875,  for  the  lossy  array  5?(e®ff)  is  continuous  (it.  drops  down 
to  -12.7,  not  shown  in  the  figure,  before  increasing  again)  with  a  rather  curious  ”V-shaped”  dip  at 
around  kd  =  0.8013.  In  the  bandgap  region  the  behavior  of  3(e®ff)  is  very  different  for  the  lossless 
and  lossy  arrays.  For  the  lossy  array,  Qf(e®ff)  is  negative  in  the  entire  bandgap  region,  decreasing  to 
-6.7  at  kd  =  0.8066  before  increasing  to  0  at  kd  =  0.8078.  In  contrast,  for  the  lossy  array  3(e*ff) 
first  decreases  to  -3.1  at  kd  =  .80,  then  jumps  suddenly  to  around  5  and  increases  rapidly  to  19  at 
kd  =  0.80875,  and  then  decreases  to  0.06  at  kd  =  0.9. 

The  effect  of  loss  on  the  behavior  of  &(/x*ff)  is  seen  from  Fig  21  to  be  significant  only  in  the 
bandgap  region,  but  the  effect  of  loss  is  much  less  pronounced  than  it  is  on  The  general 

shape  of  )  for  the  lossless  and  lossy  arrays  is  similar,  but  the  peak  value  for  the  lossless  array 
is  about  15  as  compared  with  about  6  for  the  lossy  array.  There  is  much  more  of  an  effect,  of  loss 


on  $(/irff)  than  on  For  the  lossless  array  is  positive  throughout  with  a  peak  value 

of  7.8  reached  rather  abruptly  at  kd  =  0.790  decreasing  steadily  to  0  at  kd  =  0.8078.  For  the  lossy 
array,  3(/iJ:ff)  increases  steadily  to  a  peak  value  of  4.2  at  kd  =  0.7938,  then  decreases  rapidly  to 
-2.9  at  kd  —  0.8013,  after  which  it  increases  steadily  to  positive  values  at  kd  =  0.8038. 

8.2  DIAMOND  SPHERES 

8.2.1  LINEAR  (ID)  ARRAYS  OF  DIAMOND  SPHERES 

In  this  subsubsection  we  consider  two  linear  arrays  of  diamond  spheres,  the  first  with  the  dipoles 
normal  to  the  array  axis,  and  the  second  with  the  magnetic  dipoles  parallel  to  the  array  axis.  For 
both  arrays,  er  =  5.84 ,/ir  =  1,  a/d  =  0.45.  For  the  array  with  the  dipoles  transverse  to  the  array 
axis,  the  dispersion  diagram  for  real  (3d  is  shown  in  [20,  fig.  18]  which  we  reproduce  here  as  Fig. 
22.  In  Fig.  23  we  show  the  extension  of  this  dispersion  diagram  into  the  fast-wave  region.  The 
bandgap  in  Fig.  22,  starting  at  kd  =  2.5733  close  to  the  first  magnetic  dipole  resonance  frequency 
of  kd  =  2.787  ( ka  =  1.254),  is  shown  in  Fig.  23  to  correspond  to  a  segment  with  ^ft(3d)  =  7 r  and 
a  small  non-zero  imaginary  part.  In  the  fast-wave  region,  9i(/3d)  decreases  from  the  light  line  to 
0  at  kd  =  4.10  and  then  starts  to  increase.  (We  have  truncated  the  curve  at  kd  =  4.22.)  3?(/3d) 
steadily  increases  in  magnitude  in  the  fast-wave  region,  negative  in  the  backward-wave  region  and 
switching  to  positive  when  5R(/3d)  starts  to  increase  from  0.  In  Fig.  24  we  show  the  kd-(3d  diagram 
for  the  same  array  with  two  additional  modes  shown,  both  with  high  attenuation. 

The  kd  (3d  diagram  for  a  ID  array  of  diamond  spheres  with  the  magnetic  dipoles  parallel  to  the 
array  axis  is  given  in  [20,  fig.  25]  which  we  reproduce  here  as  Fig.  25.  The  extension  of  this  diagram 
into  the  fast-wave  region  is  shown  in  Fig.  26.  We  see  that  5R(/?rf)  begins  at  the  light  line  when 
kd  =  2.407,  increases  rather  slowly  in  the  fast-wave  region  with  Q((3d)  positive  until  kd  =  2.815 
where  it  crosses  the  light  line  and  continues  to  increase  in  the  slow-wave  region  until  it  reaches  the 
value  of  7 r  close  to  the  magnetic  dipole  resonance  frequency  of  kd  =  2.787  (ka  =  1.254).  For  this 
slow- wave  region,  3(/?d)  =  0  and  this  section  of  the  curve  is  identical  to  the  curve  shown  in  Fig. 
25.  After  9? ((3d)  reaches  n  it  stays  there  for  higher  values  of  kd  with  *3((3d)  then  increasing  so  that 
the  wave  attenuates.  Note  that  for  this  kd-ftd  diagram  there  is  a  lower  cut-off  frequency  rather 
than  the  upper  cut-off  frequencies  we  have  often  encountered. 

8.2.2  2D  ARRAYS  OF  DIAMOND  SPHERES 

In  this  subsubsection  we  consider  two  2D  arrays  of  diamond  spheres,  the  first  with  the  electric 
dipoles  in  the  plane  of  the  array,  and  the  second  with  the  electric  dipoles  perpendicular  to  the 
array  plane.  As  in  the  previous  subsection,  for  both  arrays,  er  =  5.84,  \xr  —  1,  and  a/d  =  0.45.  The 
dispersion  diagram  for  a  2D  array  of  diamond  spheres  with  the  electric  dipoles  in  the  plane  of  the 
array  and  (3d  real,  is  given  in  [20,  fig.  19], [21,  fig.  10]  which  is  reproduced  here  as  Fig.  27  and  the 
corresponding  dispersion  diagram  for  complex  (3d  is  shown  in  Fig.  28.  $l(/3d)  starts  out  very  close 
to  the  light  line,  increasing  to  n  at  kd  =  2.220  (ka  =  1.000),  somewhat  lower  than  the  first  magnetic 
dipole  resonance  frequency  of  kd  =  2.787  (ka  =  1.254).  This  section  of  the  curve  is  identical  to  the 
lower  curve  of  Fig.  27.  The  bandgap  in  Fig.  27  is  shown  in  Fig.  28  to  correspond  to  a  segment 
with  9f(/?d)  =  7 r  and  3(/ 3d )  negative.  At  the  upper  end  of  the  bandgap  at  kd  =  2.507  (ka  =  1.128), 
3(/ id )  goes  to  zero  where  it  remains  until  9?(/Jd)  crosses  the  light  line.  This  section  of  the  curve 
is  identical  to  the  upper  curve  of  Fig.  27.  As  kd  continues  to  increase,  9f(/?d)  enters  the  fast -wave 
region  and  3(/3cf)  becomes  slightly  negative  (backward  wave).  After  9?(/?d)  goes  to  zero  and  then 
starts  to  increase,  3(/3cf)  switches  to  positive  values  (forward  wave).  In  Fig.  29  we  show  the  kd-(3d 
diagram  for  the  same  array  with  four  additional  modes  shown,  all  with  high  attenuation.  Figure 
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30  is  a  detail  of  Fig.  29.  Note  that  Mode  5  begins  at  the  light  line  so  that  it  h as  a  lower  cut-off 
frequency. 

The  dispersion  diagram  for  the  2D  array  of  diamond  spheres  with  the  electric  dipoles  perpen¬ 
dicular  to  the  array  plane  and  f3d  real  is  shown  in  [20,  fig.  20], [21,  fig.  11]  and  is  reproduced  here  as 
Fig.  31,  with  the  corresponding  dispersion  diagram  for  complex  (3d  shown  in  Fig.  32.  The  behavior 
of  the  dispersion  diagram  in  the  fast-wave  region  is  similar  in  essentials  to  that  in  Fig.  28,  but 
with  3?(/ 3d )  much  larger  in  magnitude  in  the  backward-wave  region.  In  Fig.  33  we  show  the  kd-(3d 
diagram  for  the  same  array  with  five  additional  modes  shown,  all  with  high  attenuation.  Figure  34 
is  a  detail  of  Fig.  33.  Like  Mode  5  in  Figs.  29  and  30,  Mode  6  begins  at  the  light  line  so  that  it 
has  a  lower  cut-off  frequency. 

The  kd  (3d  diagram  for  a  2D  array  of  diamond  spheres  with  the  magnetic  dipoles  parallel  to  the 
array  axis  and  (3d  real  is  shown  in  [20,  fig.  26], [21,  fig.  14]  and  is  reproduced  here  as  Fig.  35,  with 
the  extension  for  complex  waves  given  in  Fig.  36.  Notice  that  the  dispersion  diagram  in  Fig.  36 
does  not  extend  into  the  fast-wave  region,  unlike  the  dispersion  diagram  for  a  ID  array  of  diamond 
spheres  with  magnetic  dipoles  parallel  to  the  array  axis,  Fig.  26.  Although  the  dispersion  diagram 
for  (3d  real,  Fig.  35,  resembles  the  corresponding  ID  dispersion  diagram  of  Fig.  25  which  does  have 
an  extension  into  the  fast-wave  region,  the  two  dispersion  diagrams  are  actually  quite  different  in 
their  behavior  at  the  light  line.  In  Fig.  35,  (3d  approaches  tangency  to  the  light  line,  whereas  in 
Fig.  25  (3d  is  incident  on  the  light  line  at  a  non-zero  angle.  The  cut-off  of  the  dispersion  diagram 
of  Fig.  35  when  (id  reaches  n  is  continued  for  increasing  kd  by  a  wave  with  JR(/3d)  =  7r  and  3(/ 3d ) 
increasing  from  0.  This  continuation  ends  when  SR(/?d)  crosses  the  light  line  at  kd  =  7r.  This  kd-(id 
diagram  displays  both  lower  and  upper  cut-off  frequencies. 

8.2.3  3D  ARRAYS  OF  DIAMOND  SPHERES 

In  this  subsubsection  we  consider  3D  arrays  of  diamond  spheres  with  the  dipoles  transverse  and 
parallel  to  the  propagation  direction.  For  transverse  dipoles  the  dispersion  diagram  for  real  (3d  is 
given  in  [20,  fig.  21], [21,  fig.  12],  reproduced  here  as  Fig.  37,  and  the  dispersion  diagram  for  complex 
(3d  is  shown  in  Fig.  38.  The  lower  branch  of  the  diagram  goes  to  7 r  at  kd  =  2.024  ( ka  =  0.911) 
considerably  below  the  magnetic  dipole  resonance  frequency  of  the  spheres  at  kd  =  2.787  (ka  = 
1.254).  It  is  seen  that  the  bandgaps  in  Fig.  37  where  there  is  no  unattenuated  traveling  wave 
correspond  in  Fig.  38  to  regions  where  3(/ 3d )  ^  0  and  hence  the  wave  is  attenuated.  The  region 
3.50  <  kd  <  3.77  is  interesting  because  there  is  a  complex  wave  ($?[/?</]  >  0)  in  the  middle  of 
the  bandgap.  In  Fig.  39  shows  the  kd  (3d  diagram  for  the  same  array  with  five  additional  high 
attenuation  modes.  A  detail  of  this  figure  is  shown  in  Fig.  40. 

The  kd  (3d  curves  for  3D  arrays  of  diamond  spheres  with  the  dipoles  parallel  to  the  array  axis 
and  (3d  real  are  shown  in  [20,  fig.  31], [21,  fig.  15],  reproduced  here  as  Fig.  41,  and  the  corresponding 
curves  for  complex  (3d  are  shown  in  Fig.  42.  It  is  seen  that  for  values  of  kd  greater  than  the  values 
where  5R(/3d)  goes  to  7 r,  3((3d)  increases  steadily  with  $t((3d)  remaining  equal  to  7r.  (The  curves 
continue  upwards  along  5?(/3d)  =  it  but  we  have  truncated  them  for  clarity  of  presentation.)  Both 
curves  have  lower  cut-off  frequencies. 

8.3  SILVER  NANOSPHERE  ARRAYS 

8.3.1  LINEAR  (ID)  ARRAYS  OF  SILVER  NANOSPHERES 

In  this  subsubsection  we  consider  ID  arrays  of  silver  nanospheres.  As  discussed  in  [20,  sec.  12. 4], [21, 
sec.  6.3]  the  dispersion  diagrams  were  computed  using  the  following  Drude  model  for  relative 
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permittivity  which  agrees  quite  well  with  the  values  of  relative  permittivity  measured  by  Johnson 
and  Christy  [70]  over  the  visible  range  of  frequencies  where  lowest  order  traveling  waves  exist: 


e  u n 

€r  —  —  —  5.45  -  0.73  ,  p  .  x 
to  cj(cj  +  i'y) 


(8.1) 


with  the  plasma  frequency  ujp  =  1.72  x  1016  and  the  loss  parameter  7  =  8.35  x  1013.  To  conform  to 
the  parameters  used  in  [71],  the  radius  a  of  the  spheres  was  chosen  to  be  5  nm  and  the  spheres  were 
embedded  in  glass  with  a  dielectric  constant  equal  to  2.56  ( k  =  1.6  x  l u/c).  With  these  parameters 
the  relative  permittivity  of  the  spheres  in  glass  can  be  written  from  the  Drude  equation  in  (8.1)  as 


6 

er  — 

fglass 


2.129 


0.0234 

(A*a/1.6)[(/ra/1.6)  +  0.00139i] 


(8.2) 


Our  formulation  in  [20], [21]  assumed  real  permittivity  and  real  (3d ,  and  so  only  the  real  part  of 
(8.2)  was  used  in  the  calculations  to  obtain  the  kd-(3d  diagram.  Here  we  use  the  complex  value  of 
eT  since  (3d  is  assumed  to  be  complex. 

The  kd-(3d  curves  for  a  ID  array  of  glass-embedded  silver  nanospheres  and  electric  and  magnetic 
dipoles  perpendicular  to  the  array  axis,  using  only  the  real  part  of  the  permittivity  and  restricting 
(3d  to  be  real,  are  given  in  [20,  fig.  28], [21,  fig.  14]  which  we  reproduce  here  as  Fig.  43,  with 
the  corresponding  curves  for  glass-embedded  silver  nanospheres  using  the  complex  values  of  the 
permittivity  and  complex  (3d  are  shown  in  Fig.  44.  The  parameter  s  is  the  free-space  distance 
between  the  spheres  so  that  d  =  s  4-  2a.  Curves  are  shown  for  s  =  1  nm  (a/d  =  0.4545)  and 
.s  =  4  nm  ( a/d  =  0.3571).  The  slow-wave  region  portions  of  the  curves  for  yt((3d)  from  close  to 
the  light  line  to  close  to  7 r  are  quite  similar  for  the  real  and  complex  permittivity  cases.  Close  to 
the  light  line,  however,  the  behavior  of  the  complex  curves  is  quite  different  from  what  it  is  for  the 
real  permittivity  curves.  Whereas  (3d  for  the  real  permittivity  curves  decreases  along  the  light  line 
in  the  slow-wave  region,  in  the  complex  permittivity  dispersion  curves  3?(/?d)  crosses  the  light  line 
into  the  fast-wave  region,  decreases  to  zero,  and  then  increases  again  to  meet  the  light  line  where 
the  curves  end.  The  fast-wave  region  traveling  waves  are  significantly  attenuated.  The  propagation 
constants  also  differ  significantly  for  3?(/?d)  close  to  7r.  For  the  real  permittivity  nanospheres  there 
is  a  sharp  cut-off,  while  for  the  complex  permittivity  nanospheres  for  very  low  frequencies  there  is 
a  strongly  attenuated  wave  that  gradually  becomes  less  attenuated  as  the  frequency  increases  and 
a  region  with  small  attenuation  begins. 

The  kd  (3d  curves  for  a  ID  array  of  glass-embedded  silver  nanospheres  and  electric  dipoles 
parallel  to  the  array  axis  with  only  the  real  part  of  the  permittivity  used  and  (3d  real  are  given  in 
[20,  fig.  32], [21,  fig.  20]  which  we  reproduce  here  as  Fig.  45.  The  corresponding  curves  for  glass- 
embedded  silver  nanospheres  using  complex  permittivity  and  complex  (3d  are  shown  in  Fig.  46. 
The  addition  of  loss  and  the  extension  of  the  dispersion  diagram  into  the  fast-wave  region  results 
in  marked  changes  from  the  lossless  dispersion  diagram.  The  low-frequency  portion  of  the  lossless 
diagram  very  close  to  the  light  line  is  eliminated  and  replaced  by  a  small  “fish-hook”  section  in  the 
fast-wave  region.  This  fish-hook  region  is  shown  in  detail  in  Fig.  47.  Note  the  discontinuity  of  the 
real  and  imaginary  parts  of  (3d  at  the  light  line.  Following  this  fast-wave  fish  hook  there  is  a  long 
positive-slope  slow- wave  region  with  Q?(/3rf)  small.  Then  ^i((3d)  begins  to  increase  steadily.  First 
there  is  a  frequency  range  for  which  $i((3d)  remains  fairly  close  to  7 r,  after  which  with  increasing 
frequency  %i((3d)  decreases  as  a  negative-slope  wave  to  the  light  line  with  3f(/?d)  steadily  increasing 
to  very  large  values  (not  shown  in  the  figure). 

As  an  example  of  the  various  kd-(3d  diagrams  that  are  obtained  when  alternative  branches  of 
the  dilogarithm  and  trilogarithm  functions  are  employed,  in  Fig.  48  we  show  the  curves  of  Fig.  46 


(in  black)  along  with  a  number  of  curves  (in  color)  obtained  with  different  choices  of  the  parameters 
m2  =  0,±1  and  m3  =  0,  ±1,  (m2,  m3)  ^  (0,0),  in  (7.13a)  and  (7.136).  Note  the  curious  behavior 
of  these  alternative  branch  curves  which  end  suddenly. 

8.3.2  2D  ARRAYS  OF  SILVER  NANOSPHERES 

The  kd~(3d  curves  for  2D  arrays  of  glass-embedded  lossless  silver  nanospheres  with  electric  and 
magnetic  dipoles  perpendicular  to  the  array  axis,  the  electric  dipoles  parallel  to  the  array  plane, 
only  the  real  part  of  the  permittivity  used,  and  (id  real,  are  given  in  [20,  fig.  29], [21,  fig.  17]  which 
we  reproduce  here  as  Fig.  49.  The  corresponding  curve  for  complex  (3d  and  lossy  silver  nanospheres 
and  s  =  1  nm  ( a/d,  =  0.4545)  is  shown  in  Fig.  50.  The  behavior  of  the  lossy  dispersion  diagram 
is  different  in  an  important  respect  from  that  of  the  ID  case.  Close  to  the  light  line  there  are 
two  distinct  branches  of  the  diagram.  The  plot  of  the  upper  branch  of  5f(/3d)  crosses  the  light 
line  and  decreases  asymptotically  to  zero  as  kd  increases,  while  the  imaginary  part  of  (3d  becomes 
increasingly  negative.  For  low  frequencies  the  lowrer  branch  starts  out  very  close  to  the  light  line 
just  as  it  does  for  the  lossless  case  with  3?(/ 3d )  =  0,  but  then  as  kd  continues  to  increase  there  is  a 
small  region  in  the  slow- wave  region  (kd  from  «  0.182  to  0.216)  where  there  is  a  complex  wave  with 
positive  3(/Jd).  Close  to  5R(/?d)  =  n  the  behavior  of  the  dispersion  diagram  is  similar  to  that  of 
Fig.  28,  where  for  very  low  frequencies  there  is  a  strongly  attenuated  wave  that  gradually  becomes 
less  attenuated  as  the  frequency  increases  until  kd  ^  0.190  and  a  weakly  attenuated  region  begins. 
In  Fig.  51  we  show  the  kd-(3d  diagram  for  the  same  array  but  with  multiple  modes  displayed. 
The  corresponding  dispersion  diagrams  for  the  larger  sphere  separation  distance,  s  =  4  nm,  are 
similar  in  all  important  respects  to  those  of  the  s  =  1  dispersion  diagram  apart  from  an  expected 
displacement  of  the  curve  to  higher  frequencies,  and  are  not  shown  here. 

For  comparison  with  Fig.  50,  in  Fig.  52  we  show  the  extension  of  the  kd-(3d  diagram  of  Fig.  49 
for  complex  (3d  but  using  only  the  real  part  of  the  permittivity  instead  of  the  complex  permittivity 
used  to  obtain  the  dispersion  diagram  of  Fig.  50.  Instead  of  $l((3d)  increasing  from  the  light  line 
and  then  hooking  back  to  it  as  it  does  in  Fig.  50,  ^R((3d)  continues  to  increase  smoothly  to  7 r, 
joining  up  with  the  upper  branch  of  the  diagram  at  3?(/ Id )  =  1.09.  Also,  almost  all  the  portion  of 
the  dispersion  diagram  of  Fig.  50  down  along  5R(/?d)  «  n  is  eliminated.  The  small  branch  of  5(/?d) 
going  to  the  left  of  the  Q(/3d)  axis  at  around  kd  =  0.2  corresponds  to  the  very  small  interval  of 
yi((3d)  along  5R(/3rf)  =  7r  in  this  range  of  kd. 

The  kd  (3d  curves  for  2D  arrays  of  glass-embedded  lossless  silver  nanospheres  with  electric  and 
magnetic  dipoles  perpendicular  to  the  array  axis  and  the  electric  dipoles  perpendicular  to  the  array 
plane  are  given  in  [20,  fig.  29], [21,  fig.  29]  which  we  reproduce  here  as  Fig.  53.  The  corresponding 
curve  for  lossy  silver  nanospheres  and  .s  =  1  nm  (o,/d  =  0.4545)  is  shown  in  Fig.  54.  Again  there 
are  two  branches  of  the  lossy  dispersion  diagram.  The  behavior  of  the  upper  branch  here  is  very 
similar  to  that  of  the  dispersion  diagrams  for  the  ID  array  of  silver  nanospheres  shown  in  Fig. 
44  with  5R(/?cf)  continuing  to  decrease  from  the  light  line  in  the  fast-wave  region  to  zero  before 
increasing  to  meet  the  light  line  again  where  the  curve  ends.  The  behavior  of  the  lower  branch  of 
the  lossy  dispersion  diagram  resembles  that  of  Fig.  50,  starting  out  for  low  frequencies  very  close 
to  the  light  line  just  as  it  does  for  the  lossless  case  with  f3((3d)  close  to  0,  but  then  as  kd  continues 
to  increase  entering  a  small  region  in  the  slow-wave  region  (kd  from  ^  0.324  to  0.326)  where  there 
is  a  complex  wave  with  positive  3(/ 3d ).  The  behavior  of  the  diagram  close  to  5?(/3d)  =  7r  is  similar 
to  that  shown  in  Figs.  28  and  32.  In  Fig.  55  we  show  the  kd  (3d  diagram  for  the  same  array 
but  with  multiple  modes  displayed.  For  the  larger  separation  distance,  s  =  4  nm,  the  dispersion 
diagram  shown  in  Fig.  56  is  completely  similar  to  that  of  the  smaller  separation  diagram  in  Fig. 
54  apart  from  an  expected  shift  in  the  direction  of  higher  frequency.  The  kd-(3d  curves  for  2D 
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arrays  of  glass-embedded  silver  nanospheres  with  electric  dipoles  parallel  to  the  array  axis,  only 
the  real  part  of  the  permittivity  used,  and  (3d  real,  are  given  in  [20,  fig.  21], [21,  fig.  33]  which  we 
reproduce  here  as  Fig.  57.  The  corresponding  curves  (in  black)  for  complex  /3d  and  lossy  silver 
nanospheres  are  shown  in  Fig.  58.  Just  as  for  the  ID  case  of  a  glass-embedded  silver  nanosphere 
array  with  electric  dipoles  parallel  to  the  array  axis,  there  is  a  very  marked  difference  between  the 
real  permittivity  and  lossy  dispersion  diagrams.  For  the  lossy  case,  5R {(3d)  begins  at  the  light  line  as 
in  the  real  permittivity  case  and  increases  as  a  forward  wave  with  fairly  small  positive  3(/?d)  until 
$t((3d)  approaches  7r.  Then  3(/ 3d )  steadily  increases  and  5 R((3d)  exhibits  a  curious  behavior,  first 
remaining  close  to  n  as  the  frequency  increases,  then  decreasing  to  about  2.5,  and  then  remaining 
close  to  there  with  increasing  kd.  In  contrast  to  the  ID  case  of  a  glass-embedded  silver  nanosphere 
array  with  electric  dipoles  parallel  to  the  array  axis  (Figs.  46  and  47)  there  does  not  appear  to  be 
any  fish-hook  extension  of  the  K(/3d)  curves  into  the  fast-wave  region  close  to  the  light  line. 

In  this  same  figure,  for  the  s  —  1  nm  array  we  show  k.d-(3d  curves  obtained  when  the  Schlbmilch 
series  with  Hq  }  and  in  (3.34)  are  evaluated  using  (3.14)  and  (3.16)  using  an  alternative  branch 
of  \J Pq  —  1  rather  than  the  principal  branch.  The  red  curves  are  obtained  when  —  \//?o  -  1  is  used 
instead  of  -f  >//3q  —  1  in  both  (3.14)  and  (3.16),  and  the  green  [blue]  curves  are  obtained  when 
is  used  instead  of  *F\//?o  —  1  in  just  (3.14)  [(3.16)].  Solid  and  dotted  lines  are  used  to 
indicate  3?(/?d)  and  3(/3d),  respectively.  Of  especial  interest  is  the  way  in  which  the  red  curves 
continue  the  black  curves  obtained  using  the  principal  branch  of  \/ (3$  —  1.  It  is  also  seen  that  the 
green  curves  divide  into  two  branches  just  to  the  right  of  the  light  line,  with  one  branch  for  3f(/?d) 
traveling  up  very  close  to  the  light  line,  and  the  other  branch  for  both  sR((3d)  and  3(/?d)  joining  up 
with  the  respective  black  curves.  For  the  blue  curve,  3? {(3d)  decreases  from  0.87  to  0  as  kd  increases 
from  0  to  0.21. 

8.3.3  3D  ARRAYS  OF  SILVER  NANOSPHERES 

The  kd  (3d  curves  for  3D  arrays  of  glass-embedded  silver  nanospheres  with  electric  and  magnetic 
dipoles  perpendicular  to  the  array  axis,  only  the  real  part  of  the  permittivity  list'd,  and  (3d  real, 
are  given  in  [20,  fig.  31], [21,  fig.  19],  reproduced  here  as  Fig.  59,  and  the  corresponding  curves 
for  lossy  nanospheres  and  complex  (3d  are  shown  in  Fig.  60.  As  for  the  ID  and  2D  cases,  the 
behavior  of  the  3D  lossy  curves  differ  strongly  from  the  corresponding  real  permittivity  curves. 
For  the  complex  permittivity  case,  as  for  the  real  permittivity  case,  the  curves  begin  close  to  the 
light  line,  but  3t((3d)  instead  of  continuing  to  n  increases  only  a  little  beyond  1.0  followed  by  an 
interval  of  increasing  kd  in  which  it  crosses  the  light  line  and  decreases  to  close  to  zero.  After  this, 
with  still  further  increase  in  frequency,  8?(/3d)  increases  again,  re-crosses  the  light  line,  and  then 
remains  quite  close  to  the  light  line.  G(/?d)  is  close  to  zero  during  the  initial  increase  of  3R(/3d), 
then  increases  very  rapidly  to  around  1.0  after  which  it  decreases  with  increasing  frequency  to  close 
to  zero  where  3 ?(/?d)  begins  its  increase  from  close  to  zero,  and  remains  increasingly  close  to  zero 
thereafter. 

The  kd  (3d  curves  for  3D  arrays  of  glass-embedded  silver  nanospheres  with  the  electric  dipoles 
parallel  to  the  array  axis,  only  the  real  part  of  the  permittivity  used,  and  (Id  real,  are  given  in 
[20,  fig.  34], [21,  fig.  22],  reproduced  here  as  Fig.  61,  and  the  corresponding  curves  for  lossy 
nanospheres  and  complex  (3d  are  shown  in  Fig.  62.  The  lossy  dispersion  diagrams  differ  greatly 
from  the  corresponding  lossless  diagrams.  Instead  of  (3d  going  from  0  to  7r  in  an  extremely  narrow 
frequency  range,  here  for  very  low  frequencies  ^ft(fld)  begins  close  to  zero  with  3(/3d)  very  large.  As 
the  frequency  increases  5R(/3d)  begins  to  increase,  crosses  the  light  line,  and  continues  to  increase 
until  close  to  n.  The  portion  of  the  curves  from  a  little  before  they  cross  the  light  line  to  where 
they  reach  close  to  n  correspond  to  the  curves  for  the  lossless  case.  The  imaginary  part  of  (3d  is 
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fairly  small  in  this  portion  of  the  dispersion  curves.  As  kd  continues  to  increase  3(/?d)  increases 
rapidly.  3?(/3d)  stays  close  to  7r  for  awhile,  and  then  decreases,  somewhat  rapidly  at  first  and  then 
much  more  slowly. 
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Figure  2:  kd  (3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  fir  =  10, 
and  a/d  =  0.45;  dipoles  perpendicular  to  array  axis;  (3d  real. 


Figure  3:  kd-(3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  eT  =  / ir  =  10, 
and  a/d  =  0.45;  dipoles  perpendicular  to  array  axis;  (3d  complex. 
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Figure  4:  kd-(3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  fix  =  10, 
and  a/d  =  0.45;  dipoles  perpendicular  to  array  axis;  (3d  complex.  Multiple  modes  shown. 
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Figure  5:  kd  /3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er 
and  a/d  =  0.45;  electric  dipoles  only;  fid  real. 
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Figure  6:  kd-0d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  €r 
and  a/d  =  0.45;  electric  dipoles  only;  0d  complex. 
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Figure  7:  kd-/3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  /ir  =  20, 
and  ajd  —  0.45;  dipoles  perpendicular  to  array  axis;  /3d  real. 


Figure  8:  kd-fid  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  =  20, 

and  a/d  —  0.45;  dipoles  perpendicular  to  array  axis;  fid  complex. 
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Figure  9:  kd-(3d  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  pr 
and  a/d  =  0.45;  dipoles  perpendicular  to  array  axis;  (3d  complex.  Multiple  modes  shown. 
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Figure  10:  kd-(id  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  eT  =  fir  =  20 
and  a/d  =  0.45;  dipole  moments  parallel  to  propagation  direction;  fid  real. 


Figure  11:  kd-(id  diagram  for  a  ID  periodic  array  of  magnetodielectric  spheres  with  er  =  /rr  —  20 
and  a/d  =  0.45;  dipole  moments  parallel  to  propagation  direction;  fid  complex. 
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Figure  12:  kd-fid  diagram  for  a  2D  square- lattice  array  of  magnetodielectric  spheres  with  fr  =  //r  = 
20  and  a/d  =  .45;  electric  dipole  moments  parallel  or  perpendicular  to  the  array  plane;  /3d  real. 


Figure  13:  kd-fid  diagram  for  a  2D  square-lattice  array  of  magnetodielectric  spheres  with  er  = 
/tr  =  20  and  a/d  =  .45;  electric  dipole  moments  parallel  or  perpendicular  to  the  array  plane;  fid 
complex. 
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Figure  14:  kd-/3d  diagram  for  a  2D  square-lattice  array  of  magnetodielectric  spheres  with  er  = 
Hr  =  20  and  a/d  =  .45;  electric  dipole  moments  parallel  or  perpendicular  to  the  array  plane;  fid 
complex;  detail. 
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Figure  15:  kd-(3d  diagram  for  a  2D  square-lattice  array  of  magnetodielectric  spheres  with  er  = 
/tr  =  20  and  a/d  =  .45;  electric  dipole  moments  parallel  or  perpendicular  to  the  array  plane;  (id 
complex.  Multiple  modes  shown. 


Figure  16:  kd  (id  diagram  for  a  2D  square-lattice  array  of  magnetodielectric  spheres  with  et  — 
Hr  =  20  and  a/d  =  .45;  electric  dipole  moments  parallel  or  perpendicular  to  the  array  plane;  (id 
complex.  Multiple  modes  shown;  detail. 
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Figure  17:  kd-(3d  diagram  for  a  3D  cubic-lattice  array  of  cr  =  fiT  =  20  magnetodielectric  spheres 
(dipole  moments  normal  to  the  propagation  direction)  with  a/d  =  .45;  (3d  real. 


Figure  18:  kd-(3d  diagram  for  a  3D  cubic-lattice  array  of  fr  =  /xr  =  20  magnetodielectric  spheres 
(dipole  moments  normal  to  the  propagation  direction)  with  a/d  =  .45;  (3d  complex.  Multiple  modes 
shown. 
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Figure  19:  kd-(3d  diagram  for  3D  cubic-lattice  arrays  of  (er  =  20.0,  fir  =  20.0  4-  0.04i),  (er  - 
20.0,  fir  -  20.0  +  0.4i),  magnetodielectric  spheres  (dipole  moments  normal  to  the  propagation 
direction)  with  a/d  —  .45;  /3d  complex. 
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Figure  20:  kd-0d  curves  for  3D  cubic-lattice  arrays  of  magnetodielectric  spheres;  1)  eT  = 
(13.8,0.0),  fiT  =  (11.0,0.0);  2)  et  =  (13.8,0.1),  /zr  =  (11.0,0.0);  dipole  moments  perpendicular  to 
array  axis;  a/d  —  0.4. 


Figure  21:  Effective  permeability  and  permeability  of  3D  cubic-lattice  arrays  of  magnetodielectric 
spheres;  1)  er  =  (13.8,0.0),  Mr  =  (11.0,0.0);  2)  fr  =  (13.8,0.1),  =  (11-0,0.0);  dipole  moments 

perpendicular  to  array  axis;  a/d  =  0.4. 
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Figure  22:  kd (id  diagram  for  a  ID  periodic  array  of  diamond  spheres;  dipole  moments  normal  to 
the  propagation  direction;  a/d  =  .45;  [id  real. 


Figure  23:  led- (3d  diagram  for  a  ID  periodic  array  of  diamond  spheres;  dipole  moments  normal  to 
the  propagation  direction;  a/d  =  .45:  (id  complex. 


Figure  24:  kd (3d  diagram  for  a  ID  periodic  array  of  diamond  spheres;  dipole  moments  normal  to 
the  propagation  direction;  a/d  =  .45;  (3d  complex.  Multiple  modes  shown. 
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Figure  25:  kd-fid  diagram  for  a  ID  periodic  array  of  diamond  spheres;  magnetic  dipole  moments 
parallel  to  the  propagation  direction;  a/d  =  .45;  fid.  real. 


Figure  26:  kd-fid  diagram  for  a  ID  periodic  array  of  diamond  spheres;  magnetic  dipole  moments 
parallel  to  the  propagation  direction;  a/d  =  .45;  fid  complex. 
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Figure  27:  kd  (id  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
parallel  to  the  array  plane;  a/d  =  .45;  (id  real. 


Figure  28:  kd  -(id  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
parallel  to  the  array  plane;  a/d  =  .45;  (id  complex. 
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Figure  29:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
parallel  to  the  array  plane;  a/d  =  .45;  0d  complex.  Multiple  modes  shown. 


Figure  30:  kd-0d  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
parallel  to  the  array  plane;  a/d  =  .45;  3d  complex.  Multiple  modes  shown;  detail. 
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Figure  31:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
perpendicular  to  the  array  plane;  a/d  =  .45;  fid  real. 


Figure  32:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
perpendicular  to  the  array  plane;  a/d  =  .45;  fid  complex. 
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Figure  33:  k(l-/3d  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
perpendicular  to  the  array  plane;  a/d  —  .45;  (id  complex.  Multiple  modes  shown. 


Figure  34:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  electric  dipole  moments 
perpendicular  to  the  array  plane;  a/d  =  .45;  fid  complex.  Multiple  modes  shown;  detail. 
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Figure  35:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  magnetic  dipole 
moments  parallel  to  the  propagation  direction;  n/d  —  .45;  (3d  real. 


Figure  36:  kd-fid  diagram  for  a  2D  square-lattice  array  of  diamond  spheres;  magnetic  dipole 
moments  parallel  to  the  propagation  direction;  a/d  =  .45;  fid  complex. 
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Figure  37:  kd-fid  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  dipole  moments  normal 
to  the  propagation  direction;  a/d  —  .45;  (id  real. 


Figure  38:  kd-(3d  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  dipole  moments  normal 
to  the  propagation  direction;  a/d  =  .45;  (id  complex. 
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Figure  39:  kd-(3d  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  dipole  moments  normal 
to  the  propagation  direction;  a/d  -  .45;  fid  complex.  Multiple  modes  shown. 


Figure  40:  kd-f3d  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  dipole  moments  normal 
to  the  propagation  direction:  a/d  =  .45;  fid  complex.  Multiple  modes  shown;  detail. 
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Figure  41:  kd-(3d  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  magnetic  or  electric 
dipole  moments  parallel  to  the  propagation  direction;  a/d=  .45;  [id  real. 


Figure  42:  kd-0d  diagram  for  a  3D  cubic-lattice  array  of  diamond  spheres;  magnetic  or  electric 
dipole  moments  parallel  to  the  propagation  direction;  a/d  =  .45;  (id  complex. 


Figure  43:  kd-(3d  curves  for  a  ID  periodic  array  of  glass-embedded  silver  nanospheres;  dipole 
moments  normal  to  the  propagation  direction;  a  =  5  nm;  permittivity  real;  (id.  real. 
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Figure  44:  kd-(id  curves  for  a  ID  periodic  array  of  lossy  glass-embedded  silver  nanospheres;  dipole 
moments  normal  to  the  propagation  direction:  a  =  5  nm;  (3d  complex. 
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Figure  45:  kd-(3d  curves  for  a  ID  periodic  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm;  permittivity  real;  fid  real. 


Figure  46:  kd-(3d  curves  for  a  ID  periodic  array  of  lossy  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm;  /3d  complex. 
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Figure  47:  led- fid.  curves  for  a  ID  periodic  array  of  lossy  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm;  fid  complex;  detail. 


Figure  48:  kd-  fid  curves  for  a  ID  periodic  array  of  lossy  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm;  fid  complex;  alternative  branc  lies 
shown. 
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Figure  49:  kd  (3d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nauospheres;  electric 
dipole  moments  parallel  to  the  array  plane;  a  =  5  nm;  s  =  1, 4  nm;  permittivity  real;  (3d  real. 


Figure  50:  kd  (3d  curves  for  a  2D  square- lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  array  plane;  a  =  5  nm;  s  =  1  nm;  permittivity  complex;  >3d  complex. 
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Figure  51:  kd-(3d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  array  plane;  a  —  5  tun;  s  —  1  nm;  permittivity  complex,  (3d  complex. 
Multiple  modes  shown. 


Figure  52:  kd-(3d  curves  for  a  2D  square- lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  array  plane;  a  =  5  ran;  s  =  1  nm;  permittivity  real;  (3d  complex. 
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Figure  53:  kd0d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  perpendicular  to  the  array  plane;  a  =  5  nni;  s  =  1, 4  nm;  permittivity  real;  0d  real. 


Figure  54:  kd~0d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric- 
dipole  moments  perpendicular  to  the  array  plane;  n  =  5  nm;  s  =  1  nm;  permittivity  complex;  0d 
complex. 
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Figure  55:  kd~(3d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  perpendicular  to  the  array  plane;  a  —  5  nm;  s  =  1  nm;  permittivity  complex:  (3d 
complex.  Multiple  modes  shown. 


Figure  56:  kd  3d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  perpendicular  to  the  array  plane;  a  =  5  nm;  s  —  4  nm;  (Id  complex. 
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Figure  57:  kd-(3d  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nauospheres;  electric- 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm;  permittivity  real;  /3d  real. 


Figure  58:  kd-fid  curves  for  a  2D  square-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  —  5  nm,  /3d  complex. 
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Figure  59:  kd-fid  curves  for  a  3D  cubic-lattice  array  of  glass-embedded  silver  nanospheres;  dipole 
moments  normal  to  the  propagation  direction;  a  =  5  run;  permittivity  real;  fid  real. 
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Figure  60:  kd-/3d  curves  for  a  3D  cubic-lattice  array  of  glass-embedded  silver  nanospheres;  dipole 
moments  normal  to  the  propagation  direction;  a  —  5  nm;  permittivity  complex,  fid  complex. 
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Figure  61:  kd  (3d  curves  for  a  3D  cubic-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation;  a  =  5  nm  ;  permittivity  real;  (3d  real. 


Figure  62:  kd-(3d  curves  for  a  3D  cubic-lattice  array  of  glass-embedded  silver  nanospheres;  electric 
dipole  moments  parallel  to  the  direction  of  propagation);  a  =  5  nm;  permittivity  complex;  (3d 
complex. 
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A  COMPLEX  WAVES  SUPPORTED  BY  A 
TWO-DIMENSIONAL  SLAB 

The  key  step  of  the  analysis  we  have  performed  in  this  report  is  the  derivation  ot  the  kd-fid 
equations  for  complex  fid  by  analytically  continuing  the  corresponding  kd-fid  equations  for  real  (3d 
into  the  complex  fid  plane.  In  this  appendix  we  give  an  intuitive  justification  of  this  procedure  by 
considering  the  problem  of  complex  waves  supported  by  a  two-dimensional  magnetodielectric  slab 
above  a  ground  plane.  Although  the  structure  treated  here  is  continuous  and  not  a  discrete  periodic 
structure  as  are  the  arrays  treated  in  the  main  body  of  the  report,  the  essential  features  are  similar, 
and  the  analysis  performed  here  will  give  much  insight  into  the  analytic  continuation  procedure  used 
to  derive  the  complex  kd-fid  equation.  For  an  in-depth  treatment  of  the  subject  of  complex-plane 
representations  of  guided  waves,  see  A.  Hessel’s  chapter  on  the  general  characteristics  of  traveling 
waves,  [52,  ch.  19). 

We  consider  three  ways  of  obtaining  the  kd  fid  equation  for  the  slab  problem:  1)  solving  the 
inhomogeneous  boundary  value  problem,  2)  solving  the  homogeneous  boundary  value  problem,  and 
3)  integrating  over  the  polarization  inside  the  slab.  Throughout  the  appendix  we  assume  that  there 
is  no  variation  of  fields  and  sources  with  y  and  treat,  only  A- waves,  so  that  Ex  =  Ez  —  0.  As  stated 
earlier  in  the  report,  an  implicit  harmonic  time  dependence  of  exp(-icjt),  u  >  0,  is  assumed. 

A.l  INHOMOGENEOUS  BOUNDARY  VALUE  PROBLEM  FOR  THE  SLAB 

Consider  a  two-dimensional  slab  with  thickness  d  of  magnetodielectric  material  with  permittivity 
e  and  permeability  /i  above  a  perfectly  conducting  ground  plane  as  shown  in  Fig.  63.  The  slab  is 
illuminated  by  electric  line  sources  located  at  a  distance  L  or  greater  above  the  ground  plane.  The 
wave  equation  for  Ey(x,  z)  is  then 


V‘2Ey(x,  z )  +  k&Eyix,  z)  -  0,  d<x  <  L 

(A.la) 

V2Ey{x,  z)  +  k2Ey(x ,  z)  —  0,  0  <  X  <  d  . 

(A. lb) 

The  boundary  conditions  for  Ey(x,  z)  are 

Ey(0,z)  =  0 

(A. 2a) 

Ey(d+,Z)  =  Ey(d-,Z). 

(A. 2b) 

The  magnetic  field  H(x , 

z)  is  given  by 

H(;r,  z)  =  r— — V  x  Ey(x,  z)  y,  d<x<L 
l^/io 

(A. 3a) 

H(x,  z)  =  V  x  Ey(x,  z)  y,  0  <  x  <d 

1  ujfi 

(A. 3b) 

so  that 

1  OEJx ,  z)  .  T 

L  ■ d<i<L 

(A.4a) 

H!(x,2)=.1  0  <  X  <  d 

\u)[L  ax 

(A.4b) 

and 

dEy{d+,z)  dEy(d~,z)  , 

ax  -  ax  ■  *  “  "/w' 

(A.5) 

76 


1  Electric  line  sources  i 

\  ' 


i 

C0>M0 

L 

1 

\ 

d 

< 

e,/i 

\ 

i  K 

- > 

V  perfectly  conducting  ground  plane 


Figure  (53:  Slat)  geometry. 

Now  let  Ey{x,  z)  be  expressed  as  the  Fourier  transform 


Ey(x,z)=  [  E(x,p)e^zdP,  x  <  L 

J  —  OC 

(A. 6a) 

E(x,  p)  =  ^-  [°°  Ev (*-  A  e~i/3zdp,  x<L. 

J  — oo 

(A. 6b) 

Then  it  follows  from  (A.1)-(A.5)  that  E(x,p)  satisfies  the  following  equations: 

d2E{x{3)  +{k2  02)E{xJ3)_Oi  d<x<L 

OX 1 

(A. 7a) 

d2E(x  P)  +  (fc2  _  0l)E^  ft  _  0,  0  <x<d 
ox. 

(A. 7b) 

E(0J3)  =  0 

(A. 8a) 

E(d+,p)  =  E(d~yp) 

(A.  8b) 

and 

d  E{d+,p  d  E{d~,p 

^ r  dx  dx 

(A.9) 

From  (A. 7b)  and  (A. 8b) 

E(x,  p)  =  A(p)  sin (fc2  -  /32)1/2x,  0  <  x  <  d 

(A. 10) 

and  from  (A. 8a) 

E(x,P)  =  B{P)e^k o  “  P2)l/2A  +  C(P)  e  d(fco  P)  '  At  d<x<L 

(A. 11) 
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(A.12) 


where  A{0).  B(0),  and  C{0)  are  to  be  determined,  so  that  from  (A. 8b)  and  (A.9) 
B(/3)e‘[(^o  -  P2)'/2d]  -f  C(/3)e-*^°  -  d2)'/2d]  =  A(/?)sin[(A:2  —  02)x/2d} 


and 


^r(k2-02)^2 


B(/J)e'l(fc0  ~  P2)',2d]  _c(/3)e~l\(ka  ~  Pl)XI~d] 


=  A(/3)(fc2  -  ^2)1/2cos[(fc2  -  /?2)1/2d]  .  (A. 13) 

Asstime  that  in  (A.ll)  we  can  identify  C{0)e~x^ko  ~  P2)l/~x)  with  the  field  of  the  line  sources  in 

d  <  x  <  L  incident  on  the  slab,  and  B(0)cl^kx>  ~  P‘)XI2x\  with  the  field  scattered  from  the  slab  in 
the  region  x  >  d  so  that 


Einc(x,0)  =  C(fi)e~'l^ko- d  <  x  <  L  (A.14a) 


Esc(x,P)  =  B(/3)Mko- P2)l/2x\  x  >  d  .  (A. 14b) 

Then  (kfi  —  01)1^2  should  be  either  positive  real  or  positive  imaginary  (note  that  at  this  stage  of 
the  analysis  /?  is  real).  Let 

70  =  (*g  -  P2)1/2  (A-15a) 

and 

1  =  {k2-02)x'2  (A. 15b) 

with  7o  positive  real  or  positive  imaginary  so  that  (A.ll)  and  (A. 13)  can  be  written  as 

B(/3)  e^°d  +  C{0)  e_i^°d  =  A{0)  sin  7d  (A.16a) 


and 

B(0)  e^°d  -  C{0)  e~^°d  =  —  —  (°SV--  . 

lMr70 

(A. 16b) 

Letting 

m 

a('  )  C{0) 

(A. 17a) 

h(R,  B{0\ 
b (/?)  C{0) 

(A. 17b) 

(A.  16)  takes  the  form 

b(/3)  —  a(/3)  sin  7 d  =  — e 

(A. 18a) 

b(/3)  _  nS :  cos  7 d  = 

i/ir70 

(A. 18b) 

from  which 

b{0)ePod  +  e~[^d 

a(/j)  =  .  , 

sin  7a 

(A.19) 

with 

tan  7 d  4-  . 

b(0)  _  e-2l^d  lW0  . 

tan  7  d  —  7 - 

l^r70 

(A. 20) 
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Then  from  (A. 6a)  and  (A.  14b) 

Ev,sc(x,  z)  =  r  b(0)C(0)  e^z  d 0 

J  —  OG 
7 

tari7d-h: -  r  y  „  . 

C{0) - “  2d)l  +  Qz\  d/3,  *>d. 

-o  x  1 


-oo  tan  7^  —  7 


(A. 21) 


l/^rTO 


If  we  close  the  infinite  ft  integration  in  the  complex  ft  plane,  we  can  find  discrete  propagating  waves 
at  the  zeroes  of  the  denominator  in  (A. 21) 

7  d 


with 


and 


tan  'yd  =7  , 

l//r70« 

7d=[(Ard)2-(/3d)2]1/2 
70  d  =  [(M)2  -  m2\x'2 


(A. 22) 

(A.  23a) 
(A.23b) 

(kd)2  =  w2pe  =  (fco)Vrfr  •  (A. 23c) 

(Note  that  C(0),  because  of  its  identification  with  the  coefficient  of  the  field  of  the  line  sources 
incident  on  the  slab,  has  no  poles.)  The  zeroes  of  (A. 22),  the  kd-0d  equation,  can  be  either  real 
or  complex  depending  on  whether  they  are  located  along  the  real  axis  or  in  the  half-plane  of  the 
complex  0  plane  formed  by  closing  the  contour  of  integration  at  |/3|  =  oo.  Thus  the  propagation 
constants  of  the  complex  propagating  waves  supported  by  the  slab  are  obtained  in  effect  by  analyt¬ 
ically  continuing  the  kd-(5d  equation  for  real  fid  into  the  complex  (i  plane,  exactly  the  procedure 
we  have  employed  in  the  analysis  of  the  main  body  of  the  report.  It  is  important  to  note  that  the 
original  integration  along  the  real  0  axis  is  equal  not  only  to  the  sum  of  the  residues  at  the  poles 
of  the  kd  0d  equation  enclosed  by  closing  the  contour  of  integration,  but  also  to  the  contributions 
of  the  integrations  along  and  around  the  branch  cuts  of  the  kd—0d  equation.  Ihese  branch-cut 
integrals  yield  continuous  radiation  spectra  of  space  waves  that  can  be  excited  by  external  source's. 

A. 2  HOMOGENEOUS  BOUNDARY  VALUE  PROBLEM  FOR  THE  SLAB 

Since  the  kd-0d  equations  in  the  main  body  of  the  report  are  obtained  with  no  consideration  of 
external  sources,  we  now  solve  the  homogeneous  boundary  value  for  the  same  slab.  This  can  be 
done  simply  by  setting  C{0)  =  0  in  (A.  16)  since  C(0)  is  the  coefficient  of  the  field  of  the  line  sources 

incident  on  the  slab.  Thus  , 

B(0)  e1"1'0^  =  A{0)  sin  7 d  (A. 24a) 

aUd  BIB)  el^od  =  ( A.24b) 

Hbio 

from  which  we  immediately  obtain  the  kd- (id  equation  for  the  homogeneous  slab  boundary  value 
problem 

J-=tan7d  (A.25) 

i//r7o« 

which  is  identical  to  the  kd-0d  equation,  (A.22),  we  have  obtained  for  the  inhomogeneous  bound¬ 
ary  value  problem.  It  should  be  noted  that  although  we  have  here  obtained  the  solution  to  the 
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homogeneous  slab  boundary-  value  problem  by  specializing  the  solution  to  the  inhomogeneous  slab 
boundary  value  problem,  it  could  have  been  obtained  directly  without  reference  to  the  inhomo¬ 
geneous  slab  problem  by  treating  the  homogeneous  slab  problem  as  an  open  waveguide  problem 
and  assuming  a  propagating  mode  supported  by  the  waveguide.  The  advantage  of  proceeding  as 
we  have,  however,  is  that  the  procedure  of  obtaining  the  dispersion  equation  for  complex  0  by 
analytically  continuing  the  dispersion  equation  for  real  0  is  then  clearly  understood  as  a  simple 
consequence  of  closing  the  infinite  real  0  integration  in  the  complex  0  plane. 


A. 3  INTEGRATION  OVER  POLARIZATION 

The  derivation  of  kd-0d  equations  in  the  main  body  of  the  report  is  based  on  obtaining  the  electric 
and  magnetic  field  at  the  center  of  a  reference  sphere  by  summing  the  scattered  field  contributions 
from  all  the  other  spheres  in  the  array,  and  not  on  solving  a  boundary  value  problem.  It  is  therefore 
of  interest  to  show  that  the  kd-0d  equations  for  the  slab  that  we  have  obtained  by  solving  the  slab 
boundary  value  problems  can  also  be  obtained  by  the  procedure,  analogous  to  that  used  in  the 
array  analysis,  of  integrating  over  the  polarization  in  the  slab  to  get  the  fields  at  z  =  0  and  then 
equating  them  to  Ey(x,  0).  We  consider  the  field  of  a  single  propagating  wave  in  the  slab  which 
from  (A. 10)  and  (A. 23a)  takes  the  form 

Ey(x,  z)  =  A(0)  sin7X  e*'^2,  0  <  x  <  d,  0  real.  (A. 26) 

The  analysis  is  simpler  if  the  ground  plane  is  dropped  and  the  slab  extended  to  x  =  —  d.  For 
simplicity  we  will  assume  here  that  fir  —  1.  Then  there  is  no  magnetic  polarization  and  the  only 
source  is  the  electric  polarization 


Py(x,  z)  =  (e  —  eo)A(0)  sin7x  chiz,  -d  <  x  <  d,  -oo  <  z  <  oc. 

The  field  everywhere,  written  as  an  integral  of  the  polarization  P 7(x,  z),  is  given  by 

Ey(x,  z)  =  J  J  I  Py(x\  z')Hq  [fco\/(z-z')2  +  (*  — 2')2]  da;'d^, 

=  J°°  j'1  sin7x,ei^2'  //,!  [W(*  -  *')2  +  (*  ~  *')2]  d*'dz'  • 

At  z  =  0 


E„(x,0)  =  ifco(^  J1a(0)  J*  sm'yx'dx'  J™  cos 0z'  H&  [fc0N/(x  -  x')2  +  z/2]  dz' 
Referring  to  [58,  6.677(8) (9)] 

7o|x -  x'| 


J  cos 0z  Hi)  [fco \/(x  -  x')2  +  2,2j  dz  = 


70 


(A. 27) 


(A.28) 


(A. 29) 


(A. 30) 


where  as  usual 

^0  =  (fcj*  -  /32)1/2,  70  positive  real  or  positive  imaginary.  (A. 31) 

Substituting  (A. 30)  in  (A. 29)  and  performing  the  integration  with  respect  to  x'  we  obtain 


P„(x,  0)  =  sin7ox  (7cos7d  -  i7osin7d)  +  7osin7x 

70  L 
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(A.32) 


— : - sin  70#  (7  cos  7^  “  Ho  7 d)  4-  A(/3)  sin  7X 

7o 


But  from  (A. 26) 


(A. 33) 


Ey(x,  0)  =  A((3)  siri7x  . 


Hence,  equating  (A.32)  and  (A. 33) 


7  cos  7  d  —  i7o  sin  7  d  =  0 


(A. 34) 


or 


(A. 35) 


identical  with  the  kd-(3d  equations,  (A. 22)  and  (A. 25),  that  we  have  obtained  by  solving  the 
inhomogeneous  and  homogeneous  slab  boundary  value  problems,  respectively,  when  /*r  =  1*  Note 
that,  just  as  in  the  periodic  case  with  the  divergent  summations  for  complex  propagating  waves, 
if  we  integrate  the  polarization  of  the  complex  wave  as  in  (A. 28)  the  exponentially  increasing 
polarization  as  2  ->  +00  or  -00  depending  on  whether  3(/3)  <  0  or  >  0,  respectively,  makes  the 
integral  diverge. 
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B  DIVERGENCE  OF  COMPLEX- WAVE  FIELDS  OBTAINED 
FROM  DIRECT  SUMMATION  OF  DIPOLE  FIELDS  OF 
PERIODIC  ARRAYS 

In  this  appendix  we  address  the  issue  of  divergent  dipole-field  summations  for  complex  waves, 
namely  that  when  a  traveling  wave  supported  by  a  periodic  array  of  dipole  inclusions  has  a  complex 
propagation  constant,  the  field  incident  on  a  reference  element  of  the  array  obtained  by  summing  the 
fields  produced  by  all  the  other  array  elements  becomes  infinite.  Although  the  infinite  summation 
can  be  properly  handled  (as  we  do  in  deriving  the  kd  (3d  equation  for  complex  (3d  by  starting  with 
the  corresponding  summation  when  the  propagation  constant  is  real  and  rewriting  the  summation 
in  a  form  that  can  be  analytically  continued  into  the  complex  (3d  plane),  nevertheless  the  divergence 
of  the  summation  evaluated  at  complex  (3  remains  and  deserves  explanation.  Of  course,  a  simple 
mathematical  reason  for  the  divergence  of  the  summation  is  the  exponentially  increasing  terms 
of  the  summation  as  z  — ►  — oo  (for  complex  waves  decaying  in  the  +z  direction)  caused  by  the 
non-zero  imaginary  part  of  the  complex  (3d ,  but  we  will  show  here  that  a  deeper  physical  insight 
into  this  divergence  can  be  obtained  by  considering  not  only  a  complex  traveling  wave  itself  but 
also  the  sources  needed  to  excite  it.  While  it  may  at  first  seem  inconsistent  to  be  discussing  sources 
in  a  report  that  centers  on  obtaining  as  much  information  about  a  traveling  wave  as  is  possible 
by  treating  only  the  homogeneous  problem  without  any  explicit  consideration  of  sources,  it  will  be 
seen  that  the  introduction  of  sources  that  eventually  recede  to  an  infinite  distance  is  required  to 
elucidate  the  divergences.  In  particular,  we  want  to  explain  physically  why  the  infinite  summation 
diverges  for  complex  propagation  constants,  whereas  the  analytic  continuation  of  the  summation 
from  real  to  complex  propagation  constants  is  well  behaved. 

The  fields  of  the  traveling  waves  with  either  real  or  complex  propagation  constants  (3  satisfy 
Maxwell’s  equations.  Thus,  if  all  sources  of  these  traveling  waves  are  taken  into  account,  the 
fields  produced  by  all  the  sources  must  equal  the  fields  of  the  traveling  waves  whether  or  not  the 
propagation  constants  of  the  traveling  waves  are  real  or  complex.  Since  the  infinite  summation  of 
the  dipolar  fields  of  a  traveling  wave  produced  by  the  dipoles  of  the  array  diverges  if  the  propagation 
constant  of  the  traveling  wave  is  complex,  it  follows  that  there  must  be  sources  at  infinity  other 
than  the  dipoles  that  produce  fields  to  cancel  the  divergence  and  to  yield  the  fields  of  the  complex 
traveling  wave.  Moreover,  since  the  infinite  summation  of  the  dipolar  fields  of  a  traveling  wave 
with  real  (3  converges  to  the  fields  of  that  traveling  wave,  any  sources  at  infinity  (specifically,  at 
z  —  —  oo  assuming  that  the  power  flow  of  the  wave  is  in  the  +z  direction)  needed  to  excite  the  real 
(3  traveling  wave  contribute  negligibly  to  the  fields  of  that  traveling  wave  in  bounded  space.  The 
remainder  of  this  appendix  substantiates  this  argument  with  supporting  derivations. 

Suppose  we  want  to  excite  only  the  fields  of  a  single  traveling  wave  (decaying  in  the  -\ -z  direction 
for  complex  (3  or  with  power  flow  in  the  +z  direction  for  real  (3)  in  the  half  space  of  a  semi-infinite 
array  of  dipole  scatterers  located  at  z  =  0,  d,  2d,  3d,  •  •  •  ,  oo.  An  equivalence  principle  can  be 
derived  that  says  we  can  do  this  by  applying  surface  electric  and  magnetic  currents  on  the  plane 
2  =  0"  equal  to  -z  x  H,  and  zxEe,  respectively,  where  Er  and  H,  are  the  fields  of  the  sources  in 
the  half  space  z  <  0. 

For  ID  arrays,  the  derivation  of  this  equivalence  principle  involves  taking  the  infinite  array  of 
dipoles  carrying  the  single  traveling  wave  and  applying  the  Kottler-Franz  integral  formula  [38,  sec. 
8. 14], [72]  to  a  closed  surface  consisting  of  the  2  =  0“  infinite  plane  and  the  hemispherical  surface 
of  infinite  radius  in  the  region  z  >  0  that  caps  the  infinite  plane.  The  sources  on  this  closed  surface 
are  the  equivalent  electric  and  magnetic  surface  currents,  nxH  and  -n  x  E,  respectively,  where  n 
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is  the  outward  normal  to  the  surface.  For  proper  complex  waves,27  the  fields  of  the  traveling  waves 
decrease  exponentially  in  magnitude  as  the  radius  approaches  infinity.  Thus,  the  exponentially 
decreasing  n  x  H  and  —  n  x  E  surface  current  sources  at  infinity  do  not  contribute  to  the  fields 
inside  the  closed  surface  (as  determined  by  the  Kottler-Franz  integration)  and  only  surface  currents 
on  an  effectively  finite  extent  of  the  z  =  0_  plane  contribute  source  fields.  Moreover,  the  Kottler- 
Franz  integration  of  the  part  of  the  equivalent  surface  current  produced  by  the  dipole  sources  in  the 
half  space  2  <  0  is  identically  zero  (38,  sec.  8.14], [72]  and  thus  -z  x  H  and  z  x  E  can  be  replaced 
by  -z  x  H,  and  z  x  EP. 

For  2D  arrays,  the  same  argument  applies  using  a  semi-circular  infinite  cylindrical  surface  to 
cap  the  2  =  (T  plane.  The  axis  of  the  infinite  cylinder  is  in  the  plane  of  the  2D  array  perpendicular 
to  the  direction  of  propagation  of  the  traveling  wave,  and  the  equivalent  surface  currents  at  the 
ends  of  the  cylinder  contribute  negligible  fields  compared  to  the  fields  from  the  surface  currents  on 
the  2  =  0“  plane  as  the  length  of  the  cylinder  approaches  infinity.  For  3D  arrays  (which  have  no 
improper  waves),  semi-infinite  planes  at  x  =  ±oc  and  y  =  ±oo,  and  an  infinite  plane  at  2  =  +00 
can  be  used  to  cap  the  2  =  0“  plane.  The  contribution  to  the  fields  at  bounded  values  of  2  from  the 
Kottler-Franz  integration  on  the  infinite  2  =  +00  plane  can  be  made  to  go  to  zero  by  inserting  a 
small  imaginary  part  of  0  t.o  make  the  fields  attenuate  to  zero  as  2  — ►  +00.  1  hen  the  contribution 
to  the  fields  from  the  Kottler-Franz  integrations  on  the  semi-infinite  x  =  ±oc  and  y  =  ±00  planes 
are  also  negligible  compared  to  the  fields  from  the  surface  currents  on  the  2  =  0“  plane.  (The 
insertion  of  a  small  imaginary  part  of  0  to  attenuate  the  fields  to  zero  as  2  approaches  ±00  for  3D 
arrays  is  also  needed  to  assure  the  convergence  of  the  infinite  3D  summation  of  dipole  fields;  see 
Footnote  11  in  the  Introduction.) 

Consequently,  for  ID,  2D,  and  3D  semi-infinite  arrays,  the  fields  of  the  traveling  wave  in  the 
half  space  2  >  0  can  be  expressed  as  the  sum  Edm  X^=o  Hdn]  of  the  fields  of  the  dipoles 

plus  the  fields  [Es.  H.,]  of  the  surface  currents  in  the  plane  2  =  CT.  In  particular,  for  the  electric 
field 

E  =  ^Edn  +  E,„  2  >  0 .  (B.l) 

n=0 

Now  move  the  entire  half  space  of  dipoles  with  their  plane  of  surface  currents  to  the  left  along  the 
-2  axis  a  distance  z0  =  Nd ,  where  N  is  a  positive  integer,  while  keeping  the  coordinate  system 
fixed.  Then  the  expression  for  the  electric  field  of  the  traveling  wave  in  (B.l)  becomes 

OO 

Eo  =  E<ion  4*  E<so,  z  >  —zq  .  (B.2) 

n=-N 

27 For  outgoing  improper  waves,  that,  is,  the  “leaky  waves”  (the  waves  in  fig.  65  witli  right-pointing  arrows),  the 
exponentially  increasing  fields  at  infinity  contribute  to  the  Kottler-Franz  equivalent  surface-current  integration.  Be¬ 
cause,  however,  these  surface  currents  represent  fields  propagating  away  from  the  axis  of  the  array,  their  contribution 
to  the  fields  inside  the  closed  surface  is  no  larger  than  that  of  the  fields  produced  by  the  surface  currents  on  the 
2  =  0“  plane  and  thus  their  contribution  can  be  included  by  adding  extra  finite  surface  currents  on  the  2  =  0  plane. 
Moreover,  even  thougli  the  surface  currents  on  the  z  =  0  plane  increase  exponentially  as  the  transverse  radius  p 
approach*®  00,  they  can  be  diminished  slowly  to  zero  as  a  function  of  p  with  negligible  effect  on  the  fields  well  within 
the  transverse  radius  of  the  undiminished  surface  currents.  For  incoming  improper  waves  (the  waves  in  Fig.  65  with 
left-pointing  arrows),  the  integration  over  the  hemispherical  surface  adds  fields  with  spatial  dependence  e‘^2  + 
to  [E,.H,].  Nevertheless,  these  extra  fields  contribute  terms  that  are  finite  in  equation  (B.5)  ami  thus  do  not  change 
the  conclusions  drawn  from  (B.5). 

Because  outgoing  improper  complex  waves  (leaky  waves)  decaying  for  z  >  0  can  effectively  be  excited  by  sources  at 
2  =  0“,  such  leaky  waves  often  represent  the  dominant  fields  in  certain  regions  of  space  in  the  solutions  to  waveguide 
excitation  problems.  In  contrast,  the  incoming  improper  waves  require  excitation  surface  currents  on  both  the  2  =  0 
plane  and  its  capping  surface.  Thus,  incoming  improper  waves  do  not  generally  represent  a  significant  contribution 
to  the  fields  in  the  solution  to  waveguide  excitation  problems  [64]. 
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For  z0  (or  equivalently  N)  large  enough,  the  electric  field  of  the  surface  currents  (now  on  the 
2  —  ( — 2q)  —  plane)  at  a  large  distance  from  the  plane  will  have  the  general  far-field  behavior 
Es0  «  F.,eifc2° /f{zo)  where  /(z0)  =  zo,  and  1  for  ID,  2D,  and  3D  arrays,  respectively.  Fs 
is  not  a  function  of  zq  and  the  approximation  becomes  exact  as  z0  — >  oo.  Thus,  as  zq  gets  large, 
equation  (B.2)  can  be  rewritten  as 


Eo  =  EjOn  + 
n=-N 


Z  >  Zq  , 


(B.3) 


However,  as  zq  — ►  oo  the  fields  [Eo,  Ho]  of  a  traveling  wave  with  a  complex  propagation  constant 
0  will  become  exponentially  small  at  any  finite  value  of  2.  Fortunately,  the  original  values  of  the 
fields  of  the  traveling  wave  that  existed  in  the  initial  half  space  z  >  0  can  be  restored  by  simply 
multiplying  (B.3)  by  the  exponential  e~'^z°  to  get 


„-i/3zo  ., 

e  --  jkzo 


E  =  VEdnf  J  Fselkz°,  z>-z0 
^  f(z  o) 


n=-N 


(B.4) 


since  E  =  Eoe-^^  and  E(in  =  Ed0„e_i^2°.  With  0  complex  and  Q(0)  >  0  for  complex  waves 
that  decay  in  the  +z  direction,  the  factor  / f{z0)  -*  oo  as  z0  ->  oo,  and  equation  (B.4) 

becomes 

oo 

E  =  ^Edn  +  (oo)Fsel/c2°,  z  >  oo .  (B.5) 

n=— OO 

This  equation  shows  that  the  surface-current  sources  required  to  maintain  a  non-zero  complex 
traveling-wave  electric  field,  E,  at  finite  values  of  z,  produce  an  electric  field  (last  term  in  (B.5)) 
that  diverges  as  these  surface-current  sources  recede  to  z  =  — oo  (20  — >  +00).  Since  the  traveling- 
wave  electric  field  E  in  (B.5)  is  finite  for  all  finite  2,  it  follows  that  the  summation  in  (B.5)  of  the 
fields  produced  by  the  infinite  array  of  dipoles  supporting  the  complex  traveling  wave  must  also 
diverge.  Although  this  divergence  is  mathematically  obvious  because  of  the  terms  in  the  summation 
that  grow  exponentially  large  as  n  — *  -00  (as  mentioned  at  the  beginning  of  this  appendix),  we 
have  now  shown  physically  that  this  divergence  is  a  consequence  of  the  infinite  fields  produced  by 
the  surface-current  sources  at  2  =  —00  that  are  required  to  excite  a  complex  wave  that  has  finite 

fields  at  bounded  values  of  z.  „ 

If,  however,  the  traveling  wave  has  a  purely  real  propagation  constant,  the  factor  e  '•  2°/ }{zq)  — ► 
0  as  zq  — ►  00  for  ID  and  2D  arrays,  and  (B.4)  reduces  to 

OO 

E  =  ^Ed„,  2  >-oc  (B.6) 

n=— OO 

which  implies  that  the  fields  of  a  traveling  wave  with  a  real  propagation  constant  are  equal  to  the 
fields  produced  by  the  dipoles  alone,  and  that  the  fields  of  the  sources  at  2  =  -00  required  to  excite 
a  real  0  traveling  wave  do  not  contribute  to  the  fields  at  bounded  values  of  2  (because  they  are 
proportional  to  l//(zo)  — '  0  as  20  — *  00  for  ID  and  2D  arrays). 

For  traveling  waves  with  purely  real  propagation  constants  on  3D  arrays,  f(z0)  —  1  and  (B.4) 

becomes 

E  =  ^Edn  +  e1(fc_(3)ZoFs,  2  > -00  (B.7) 
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which  confirms  that  the  infinite  summation  of  unattenuated  dipole  fields  of  a  3D  array  converges 
except  for  an  oscillatory  contribution  given  by  the  last  term  in  (B.7).  This  oscillatory  term  can  be 
suppressed  in  the  infinite  summation  of  dipole  fields  by  giving  the  dipoles  a  small  attenuation  as 
2  — »  ±oo;  see  Footnote  11  in  the  Introduction. 
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C  PROPERTIES  OF  COMPLEX  WAVES  ON  UNIFORM  OR 
PERIODIC  WAVEGUIDES 

This  appendix  is  divided  into  two  sections.  In  Section  C.l,  we  summarize  the  general  properties  of 
complex  waves  supported  by  uniform  or  periodic  waveguides.  In  Section  C.2,  we  discuss  some  ol 
the  distinguishing  features  of  fast  and  slow  waves.28 

C.l  INTRODUCTION  TO  COMPLEX  WAVES 

Much  of  the  information  on  complex  waves  summarized  in  this  appendix  can  be  found  in  greater 
detail  in  the  references  [64]  and  [36].  We  shall  use  the  term  “complex  wave”  to  refer  to  a  time- 
harmonic  ,  oj  >  0)  discrete  wave  with  a  longitudinal  (z  axis)  propagation  constant  0  that  is 

neither  purely  real  nor  purely  imaginary.  The  discrete  waves  with  purely  real  propagation  constants 
(3  on  one-dimensional  and  two-dimensional  lossless,  open  waveguides  are  commonly  referred  to 
as  surface  waves,  and  waves  with  purely  imaginary  propagation  constants  0  (on  open  or  closed 
waveguides)  are  commonly  referred  to  as  evanescent  waves. 

We  begin  by  considering  uniform  or  periodic,  open  (to  free  space)  or  closed  infinite  wave¬ 
guides  composed  of  lossy  or  lossless,  reciprocal  material.  The  waveguides  can  be  one  dimensional 
(ID)  (uniform  or  periodic  in  the  2-propagation  direction  but.  not  in  the  x  and  y  directions),  two 
dimensional  (2D)  (uniform  or  periodic  in  the  z-propagation  direction  and  uniform  or  periodic  in 
either  the  x  or  y  direction  but  not  in  both  the  x  and  y  directions),  or  three  dimensional  (3D) 
(uniform  or  periodic  in  the  z-propagation  direction,  uniform  or  periodic  in  the  x  direction,  and 
uniform  or  periodic  in  the  y  direction).  Fully  periodic  ID,  2D,  or  3D  waveguides  are  formed  by  ID, 
2D,  or  3D  infinite  rectangular-lattice  arrays  of  identical  elements.  All  3D  waveguides  (uniform,  and 
periodic)  can  be  considered  closed  waveguides  with  respect  to  waves  that  are  uniform  or  uniformly 
periodic  in  the  two  rectangular  directions  transverse  to  the  z-propagation  direction  because  there  is 
no  radiation  loss  transverse  to  the  z-propagation  (longitudinal)  direction.  For  waveguides  that  are 
periodic  in  the  z  direction,  0  refers  to  the  primary  Floquet-mode  propagation  constant  satisfying 
|»(/3d)|  <  7T.  The  complete  spectrum29  of  waves,  referred  to  as  modes,  on  closed  waveguides  at  any 
one  frequency  consists  of  an  infinite  number  of  discrete  modes.  The  complete  spectrum  of  waves 
(modes)  on  lossless  uniform  (in  the  z  direction)30  open  waveguides  at  any  one  frequency  consists 
of  a  denumerable  (finite  or  infinite  depending  on  the  waveguide)  number  of  discrete  modes  plus  a 
continuous  spectrum  of  “radiation  modes.”  The  discrete  spectrum  of  modes  on  open  waveguides  is 
composed  of  all  independent  homogeneous  solutions  to  Maxwell  s  equations  that  have  exponentially 
decreasing  fields  in  the  direction  transverse  to  the  z-propagation  axis  (longitudinal  direction)  of  the 
waveguide  (and  thus  carry  finite  power  along  the  z  direction).  These  discrete  modes  are  called 
proper  waves.  The  corresponding  discrete  waves  on  lossless  open  waveguides  that  obey  Maxwell's 

28  A  good  source  for  early  references  on  the  subject  of  complex  waves  on  isotropic  structures  plus  a  treatment  that 
in  some  ways  complements  ours  is  given  in  [73].  Also  see  the  treatment  of  traveling  waves  by  A.  Hessel  in  [52,  ch. 
10],  and  the  treatment  of  leaky  waves  by  T.  Tamir  in  [52,  ch.  11]  and  by  D.  Jackson  and  A.  Oliner  in  [37,  sec.  IV]. 

29The  word  “complete”  is  used  here  in  its  mathematical  sense,  as  in  a  “complete  set  of  eigenfunctions,”  to  refer  to 
the  full  set  of  independent  solutions  satisfying  wave  equations  and  boundary  conditions  including  boundary  conditions 
of  finite  fields  at  infinity. 

30Open  waveguides  periodic  in  the  2-propagation  direction  always  have  an  infinite  number  of  discrete  higher  order 
Floquet  modes  associated  with  each  primary  Floquet  mode.  Also,  as  indicated,  there  can  be  an  infinite  number  ol 
discrete  complex  modes  (proper  waves)  on  open  waveguides  that  are  uniform  in  the  2-propagation  direction  [74],  and 
we  cannot  rule  out  the  possibility  of  an  infinite  number  of  primary  complex  Floquet  modes  on  open  waveguides  that 
are  periodic  in  the  2  propagation  direction.  There  can  also  exist  an  infinite  number  of  discrete  improper  complex 
waves  on  some  waveguides  [74]. 
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homogeneous  equations,  but  have  exponentially  increasing  fields  in  the  transverse  direction,  carry 
infinite  power  along  the  2  direction  and  thus  are  not  part  of  the  complete  spectrum  of  modes.  They 
are  called  non-spectral  or  improper  waves.  " 

It  was  proven  in  [36]  that  every  complex  (proper  or  improper)  wave  on  a  lossless,  reciprocal, 
open  or  closed,  uniform  or  periodic  waveguide  comes  as  one  of  a  quadruplet  of  complex  waves  with 
propagation  constants  0,  -0,  0 *,  and  -0*.  Two  of  these  complex  waves  decay  exponentially  in  the 
+z  direction  and  two  decay  exponentially  in  the  -z  direction.32  Because  the  complex  waves  that 
increase  exponentially  in  the  z  direction  away  from  the  sources  (assumed  finite  in  magnitude  and 
extent)  would  dissipate  an  infinite  amount  of  power  if  a  small  loss  were  inserted  into  the  material 
of  the  waveguide,  only  the  two  complex  waves  decaying  away  from  finite-power  sources  may  exist 
to  the  right  and  to  the  left,  of  the  sources.  That  is,  any  allowable  complex  waves  must  attenuate  in 
the  2  direction  away  from  the  sources.33  Assuming  0  is  the  propagation  constant  of  a  complex  wave 
decaying  in  the  +z  direction  to  the  right  of  the  sources  exciting  a  lossless,  reciprocal  waveguide,  we 
can  confine  our  attention  to  the  two  complex  waves  with  propagation  constants  0  and  - 0 *  because 
the  complex  waves  with  propagation  constants  -0,  and  0"  merely  propagate  with  decay  in  the  -z 
direction  to  the  left  of  the  sources. 

In  [36]  it  was  also  proven  that  the  electric  and  magnetic  fields  (E,  H)  of  a  complex  wave  with 
propagation  constant  0  on  a  lossless,  reciprocal  waveguide  simply  change  to  (E*,  -H*)  for  the  cor¬ 
responding  complex  wave  with  propagation  constant  —0*.  Therefore,  on  open,  lossless,  reciprocal 
waveguides,  the  two  complex  waves  decaying  in  the  +z  direction  with  propagation  constants  0  and 
-0*  are  either  both  proper  waves  (both  part  of  the  complete  spectrum  of  discrete  modes,  which 
have  fields  that  decrease  exponent  ially  transverse  to  the  2  axis  of  the  waveguide)  or  both  improper 
waves  (with  fields  that  increase  exponentially  transverse  to  the  z  axis  of  the  waveguide  and  are 
not  part  of  the  complete  spectrum  of  discrete  modes).  The  exponential  variation  of  the  far  fields 
of  these  two  complex  waves  decaying  to  the  right  of  the  sources  in  the  +z  direction  can  be  found 


directly  from  Maxwell's  equations  as 

eiKPei0z  and  e-‘«V  e-^*z 

(C.l) 

with 

0  =  5?(/?)  +  i 3(/3) ,  3(/J)  >  0 

(C.2) 

and 

k  =  iy/k2  —  /J2  ,  k  =  lj/c 

(C.3) 

31  For  evanescent  waves  (3?(/3)  =  0),  the  proper  and  improper  criteria  of  exponentially  decreasing  and  increasing 
waves  in  the  transverse  direction  are  replaced  by  outgoing  and  incoming  waves  in  the  transverse  direction.  The  limit 
case  of  a  TEM  surface  wave  (3(/i)  =  0)  with  |SJ(/3)|  =  k  is  not  a  proper  wave  because  it  carries  an  infinite  amount  of 
power  in  the  longitudinal  (2)  direction. 

32Since  the  waveguides  are  lossless,  the  exponential  decay  is  consistent  with  energy  conservation  only  if  the  total 
average  (real  part  of  the  complex  Poynting  vector)  power  flow  is  zero  for  each  proper  complex  wave  across  a  transverse 
plane  (z  =  constant  plane).  However,  a  pair  of  proper  complex  waves  with  2-propagation  constants  0  and  -0*  do 
not  generally  display  complex  Poynting  vector  orthogonality  and  can  transfer  net  complex  power  in  the  2  direction 
[75].  Still,  it  can  be  shown  that  both  E  x  H  and  real  power  flow  for  0  and  -0m  proper  complex  modes  exhibit 
orthogonality  [73].  But  unlike  opposite  decaying  identical  evanescent  waves  [76,  p.  65],  equation  (99)  in  [75]  can  be 
used  to  show  that  opposite  traveling  identical  proper  complex  waves  (with  2-propagation  constants  0  and  -0)  exhibit 
complex  power  orthogonality  and  transfer  no  net  complex  power  [73].  (As  mentioned  in  the  previous  paragraph,  the 
total  power  flow  in  the  2  direction  for  an  improper  complex  wave  on  an  open  waveguide  is  infinite  —  a  property  that 
disallows  them  from  the  complete  spectrum  of  open  waveguides.) 

33This,  however,  does  not  mean  that  a  complex  wave  is  a  proper  wave  if  it  attenuates  away  from  the  sources  in 
the  2  direction  because  improper  complex  waves  on  open  waveguides  can  attenuate  away  from  the  sources  but  are 
disallowed  as  part  of  the  complete  spectrum  because  they  carry  infinite  power  in  the  2  direction  as  a  result  of  their 
exponential  increase  transverse  to  the  2  direction. 
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where  k  is  the  transverse  propagation  constant  of  the  primary  Floquet  mode  (assuming  its  fields 
are  not  zero)  with  corresponding  2-propagation  constant  0 ,  and  p  refers  to  a  direction  transverse  to 
the  2  axis.34  If  the  sign  of  k  in  (C.3)  is  chosen  to  make  the  3(k)  >  0,  then  the  two  complex  waves 
in  (C.l)  on  the  open,  lossless,  reciprocal  waveguide  are  proper  waves.  (For  2D  open  waveguides,  we 
assume  that  the  sign  of  k  is  the  same  on  both  sides  of  the  waveguide.  For  the  more  general  case, 


see  [77].)  The  two  complex 

waves  in  (C.l)  can  also  be  written  as 

eik  '  r  and  e"ik*  '  r 

(C.4) 

with 

k  =  up  +  0i ,  r  =  pp  +  2Z 

(C.5) 

and  thus 

k  k  =  k*  k*  =  k2 . 

(C.6) 

Equation  (C.l)  or  (C.4)  shows  that  the  0  and  -0*  waves  correspond  also  to  k  and  -k*  waves.  Thus 
t  he  only  difference  between  the  full  propagation  behavior  of  the  0  and  —0*  waves  is  the  sign  of  the 
phase.  From  the  equations  in  (C.6),  we  find  that  5?(k)  -9;(k)  =  0.  Therefore,  if  one  moves  along  the 
direction  given  by  3?(k),  specifically,  r  =  r0  +  o'ft(k)  with  a  a  real  number  and  r0  a  fixed  reference 
position,  there  is  no  change  in  the  magnitude  of  the  waves  in  (C.4)  or  (C.l).  In  other  words,  S?(k)  is 
orthogonal  to  the  equi-phase  planes.  Moreover,  (C.4)  shows  that  the  magnitude  of  the  wave  varies 
exponentially  in  these  equi-phase  planes.  That  is,  in  the  far-field,  the  local  equi-phase  and  equi- 
magnitude  lines  are  orthogonal.  Consequently,  with  just  constant-magnitude  rays,  we  can  depict 
the  two  proper  complex  waves  determined  by  (C.4)  or  (C.l)  as  shown  in  Fig.  64,  and  the  two 
improper  complex  waves  determined  by  (C.4)  or  (C.l)  as  shown  in  Fig.  65.  I  he  twoheaded  arrows 
in  each  of  the  Figs.  64  and  65  along  the  rays  are  the  respective  directions  of  increasing  phase  of  the 
two  complex  waves  [the  direction  of  31(k)].  A  larger  distance  between  the  solid  constant -magnitude 
rays  indicates  a  smaller  magnitude  of  the  complex  waves.  In  the  directions  normal  to  the  rays 
(arrows),  the  phase  is  constant  and  the  fields  vary  exponentially.  This  exponential  variation  is  in 
the  direction  of  3(k)  (along  the  dotted  lines  in  Figs.  64  and  65). 

Since  the  improper  complex  waves  in  Fig.  65  are  not  part  of  the  complete  spectrum  of  modes 
on  the  open,  lossless,  reciprocal  waveguide,  one  could  ask  why  we  include  them  in  our  discussion. 
For  one  thing,  these  improper  waves  can  convert  to  proper  waves  if  the  open  waveguide  is  converted 
to  a  closed  waveguide,  for  example,  by  surrounding  it  with  a  perfectly  conducting  cylinder  [78]. 
Secondly,  these  improper  waves  still  satisfy  Maxwell’s  homogeneous  equations  everywhere  (r  <  oo) 
even  though  they  are  not  finite  as  p  -»  oo  and,  thus,  they  are  of  theoretical  interest.  More 
importantly,  however,  in  certain  limited  regions  in  the  vicinity  of  the  material-air  boundaries  of  the 
waveguide,  the  improper  wave  with  outward  propagating  energy,  depicted  by  the  arrows  pointing 
toward  the  right  (left)  in  Fig.  65  for  forward  (backward)  improper  complex  waves,  can  represent 


34  If  the  primary  Floquet  mode  decays  exponentially  in  the  p  direction,  all  the  Floquet  modes  decay  exponentially 
in  the  p  direction. 
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Figure  64:  The  two  proper  complex  waves  for 
2  >  0  (to  the  right  of  the  sources)  on  lossless, 
reciprocal,  open  waveguides. 


Figure  65:  The  two  improper  complex  waves 
for  2  >  0  (to  the  right  of  the  sources)  on  loss¬ 
less,  reciprocal,  open  waveguides. 


the  dominant  portion  of  the  fields  [64],  for  example,  on  open  dielectric-slab  waveguides  with  per¬ 
mittivity  greater  than  that  of  free  space  [74].  This  outward-energy-propagating  improper  wave 
on  an  open  waveguide  is  commonly  referred  to  as  the  improper  leaky  wave.  The  inward-energy- 
propagating  improper  wave,  depicted  by  the  arrows  pointing  toward  the  left  (right)  in  Fig.  65 
for  forward  (backward)  improper  complex  waves,  has  apparently  never  been  found  to  represent  a 
significant  portion  of  the  fields  anywhere  [64],  presumably  because  it  represents  energy  coming  from 
sources  at  p  =  oo  and  2  to  the  right  of  the  waveguide  sources;  see  Footnote  27.  Proper  complex 
waves  are  not  found  on  open  dielectric-slab  waveguides  with  permittivity  greater  than  that  of  free 
space,  but  are  found,  for  example,  on  slab  waveguides  made  with  double  negative  materials  [74]. 
Here  in  this  report  we  find  complex  waves  well  within  the  slow-wave  region  where  we  expect  they 
are  proper  waves  (see,  for  example,  Fig.  18). 

In  Section  8  of  this  report,  we  compute  numerical  values  for  the  complex  propagation  constants 
{(3d)  versus  frequency  ( kd )  on  ID,  2D,  and  3D  periodic  arrays  of  lossy  and  lossless,  reciprocal 
magnetodielectric  spheres.  For  the  lossless  open  arrays  (that  is,  for  the  ID  and  2D  arrays),  these 
complex  values  of  (3d  alone  do  not  uniquely  determine  the  proper  or  improper  nature  of  their 
associated  complex  waves.  For  example,  suppose  we  find  that  9t(/3)  >  0  and  S(/3)  >  0  at  a 
particular  frequency.  The  sign  ambiguity  in  (C.3)  allows  this  wave  to  be  either  an  improper  leaky 
wave  depicted  in  Fig.  65  with  arrows  pointing  toward  the  right,  or  a  proper  complex  wave  depicted 
in  Fig.  64  with  arrows  also  pointing  toward  the  right.  One,  of  course,  could  distinguish  between 
these  two  possibilities  with  additional  information  such  as  whether  the  waves  are  exponentially 
decreasing  (proper  wave)  or  increasing  (improper  wave)  with  transverse  distance  from  the  array; 
see  last  paragraph  of  Section  7. 

On  lossy  waveguides,  the  complex  conjugate  waves  become  two  separate  waves  with  propagation 
constants  not  equal  to  the  complex  conjugate  of  each  other  [79].  Therefore,  to  find  backward¬ 
traveling  waves  on  lossy  waveguides,  both  positive  and  negative  values  of  ^{(3)  must  be  considered  to 
include  all  the  +z  decaying  complex  waves  if%t((3)  is  restricted  to  positive  values  in  the  computations. 

C.2  FAST  AND  SLOW  WAVES 

Whether  the  discrete  wave  is  on  a  lossy  or  lossless,  reciprocal  or  nonreciprocal,  open  waveguide, 
the  equation  for  n  in  (C.3)  describes  the  transverse  far  field  of  the  wave  given  by  the  exponential 
function  elKPex@z  in  (C.l).  With  the  real  and  imaginary  parts  of  the  2-propagation  constant  written 
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as  /?  =  ft  +  i/3,,  the  transverse  propagation  constant  k  can  be  expressed  from  (C.3)  as 

k  =  i^/fc2  +  ft  -  /32  -  2iftft  (C.7) 

wliere  we  are  concentrating  on  waves  to  the  right  side  ( +z  side)  of  the  sources  so  that  ft  >  0.  The 
discrete  wave  is  usually  classified  as  a  fast  wave  it  |/?r|  <  k  and  as  a  slow  wave  it  |ft-|  >  k. 

For  ftft  ^  0,  we  see  from  (C.7)  that  k  has  a  non-zero  imaginary  part  and  thus  the  sign  in  (C.7) 
can  be  chosen  to  yield  a  proper  (or  an  improper)  complex  wave.  Note,  however,  that  we  arrive  at 
a  contradiction  if  we  assume  that  ft  =  0  and  |ft|  <  k,  because  then  k  would  be  real  and  represent 
radiation  outgoing  to  p  =  +00  (for  the  +  sign  in  (C.7))  or  radiation  incoming  from  p  =  +00  (for  the 
-  sign  in  (C.7)).  This  radiation  to  or  from  p  =  +00  would  not  be  compatible  with  an  unattenuated 
wave  (ft  =  0)  in  the  z  direction.  Therefore,  for  the  fast-wave  region  |ft|  <  k,  the  wave  must  be 
complex  (or  evanescent  if  ft  =  0).3,> 

For  ft  =  0  on  lossless,  open  waveguides  and  |ft|  >  k,  the  discrete  wave  is  a  proper  (exponentially 
decreasing  with  p)  or  improper  (exponentially  increasing  with  p)  surface  wave  depending  on  the 
sign  of  the  square  root  in  (C.7).  For  ft  >  0,  there  is  no  reason  that  a  discrete  wave  cannot  be 
complex  in  the  slow-wave  region,  that  is,  for  |ft|  >  k,  even  for  lossless,  open  waveguides.  For 
example,  if  ft  =  k  on  a  lossless  waveguide  in  (C.7),  then  k  =  0  (the  limit  case  of  a  surface  wave) 
for  ft  =  0,  but  k  is  complex  for  ft  >  0.  In  other  words,  it  is  possible  for  a  non-zero  value  of  ft 
in  a  slow  traveling  wave  (ft  >  k)  to  produce  radiation  in  the  transverse  direction  compatible  with 
energy  conservation.  This  result  has  been  confirmed  by  the  asymptotic  and  numerical  evaluations 
of  k-(3  curves  in  [64]  and  [74].  The  slow-wave  range  k  <  |ft|  <k  +  Ak(Ak>  0)  where  the  discrete 
wave  is  complex  on  a  lossless,  open  waveguide  is  said  by  Tamir  and  Oliner  [64]  to  be  the  region 
where  there  is  a  “winding  down  of  the  leaky-wave  solution.”  Beyond  this  range  (|ft|  >  k  +  Ak), 
the  wave  becomes  a  surface  wave  with  ft  =  0  on  lossless,  open  waveguides. 


“An  exception  to  this  rule  would  be  a  ID  or  2D  lossless  periodic  array  with  a  null  in  the  embedded  element 
pattern  such  that  the  fast  wave  having  a  phase  progression  corresponding  to  this  null  direction  would  see  a  radiat  ion 
“blind  spot"  [25].  Such  a  radiationless  fast  wave  would  be  a  surface  wave  having  no  primary  Floquet  mode  (for  which 
k  =  ±y/k2  -  /3? )  but  its  dominant  non-zero  Floquet  mode  still  described  in  the  far  field  by  e'Kf“t^e^  “  with  0,  =  0, 
\0r\  <  k,  and  tcfMt  imaginary  [80]. 
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